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Chapter 1 
Introduction 



The theory of the strong interaction is known to be Quantum Chromo Dynamics (QCD). It 
has all the features which are necessary for a successful description of the strong interac- 
tion. There are very important reasons why it is necessary to invest a lot of effort to solve 
QCD: 

• Validate or invalidate QCD by comparing its predictions with experiments. Results 
in the high energy regime show very good agreement with the experiments, however 
there are still many white areas with no results at all, among them is the missing 
connection between nuclear physics and QCD. 

• Validate or invalidate the Standard Model of particle physics. Even if QCD is the 
proper theory of strong interactions, it can happen that the weak and electromag- 
netic interactions are not correctly described by Standard Model. Examining weak 
decays can only be done by taking into account low-energy strong interaction effects. 
The success of the (in)validation is now mostly depends on the precision of QCD cal- 
culations. 

• Unfold the phase diagram and properties of QCD at finite temperature and baryon 
densities. In parallel with the theoretical developments intensive experimental work 
is done (and will be done) to produce and investigate the high temperature phase of 
QCD: the quark-gluon plasma. Among these investigations the major goal is to find 
signals of a first-order or second-order transition. 

Solving the above problems is known to be extremely difficult. Currently available 
methods (e.g.. weak coupling perturbation, 1/N C expansion^, string theory methods, lat- 
tice) are not able to provide us with rigorous solutions. However some of these methods 
are believed to give us very good approximations of these solutions. 

^N c is the number of colors, in QCD N c = 3. 
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1.1. Dynamical f ermions in lattice QCD 



Today the lattice technique is the only which is (or will soon be) able to calculate masses 
of hadrons, properties of low energy scattering processes, bulk and spectral properties of 
the quark-gluon plasma and many more based only on the Lagrangian of QCD. It contains 
systematical errors, however these can be quantified, therefore can be kept under control. 
In the following short introduction to the lattice technique we will highlight the role of 
"dynamical fermions" in lattice QCD. 

1.1 Dynamical fermions in lattice QCD 

Lattice QCD discretizes^l the path integral (Z) 

Z= [[dU][dxl> f ][d$ f ]exp(-S sau& -£VfD[U]rp f ) (1.1) 
J f 

on a four dimensional Euclidean lattice. The Euclidean space formalism is useful for 
obtaining the spectrum of the theory or for doing finite temperature calculations. The 
Minkowski space approach, which is necessary to investigate real time processes, is not 
available (however see flU). In Eq. Il.ll we have an integral over the gauge (U) and flavored 
fermion fields H (ipf, tpf). The S gau ge is the gauge action, the fermion action is bilinear in 
the fermion fields, so the fermion integral can be easily carried out. We end up with the 
determinant of the Dirac-operator (D[U]) under the path integral: 

Z = J [dU] det D[U] n f exp(-S gaug e), (1.2) 

with n t being the number of fermion flavors. The minimal distance on the lattice is called 
lattice spacing (a). The final results are obtained by sending the lattice spacing to zero 
together with doing some necessary renormalization. 

There are two important observations: firstly the Euclidean path integral is equiva- 
lent with a Boltzmann sum of a statistical mechanical system and secondly the discretized 
path integral in a finite volume can be put on a computer. These properties made lattice 
QCD a multidisciplinary science: it is a mixture of quantum field theory, statistical physics, 
numerical analysis and computer science. The development in computer algorithms and 
the exponential rise of the available computational capacity made lattice QCD from a toy 
model to a powerful predictive tool, giving us high precision pre- and postdictions in a 
huge number of areas. 

Nowadays the calculations are reaching the % level precision thanks to the gradual 
elimination of the so called quenching effects. Quenching means approximating the fermion 

2 Introductory materials covering the extended literature are |TJ[2j[3]|. Annual review of the field can be 
found in the lattice conference proceedings. To avoid the proliferation of citations in the introduction we 
refer to these, and cite articles only in special cases. 

3 For simplicity we have same quark masses for the different flavors in this introductory section. 
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1.1. Dynamical fermions in lattice QCD 



determinant with a constant, If independent value in Eq. 11.11 and keeping fermions only 
in the correlation functions. It is used to decrease the computational requirements, since 
taking into account the fermion determinant in the path integral (in other words dealing 
with the fermions dynamically) is a hard task. One can look on the system of ll.ll as a gauge 
system but with a highly nonlocaj^ effective action: 

S e ff = S gauge ~ «/ log det D. 

Developing efficient algorithms for such systems is nontrivial, however there is a consid- 
erable progress in the last years. 

There is a huge arbitrariness in choosing the type of the discretization, only a few re- 
quirements are to be fulfilled: eg. it should have appropriate symmetries or there should 
exist an equivalent local formulation. Universality, a well-known concept from statistical 
physics ensures that in the zero lattice spacing limit the results will not depend on the 
choice of the discretization. Since the cost of algorithms usually goes with an enormous 
power of the inverse lattice spacing, in practice it is desirable to improve the lattice ac- 
tions, that is to reduce their lattice artefacts to make the continuum extrapolations easier 
from the available lattice spacings. However one should be careful with the improvement: 
overimproving can lead to several practical problems (loss of locality, unitarity, irregular 
continuum limit, slowing down of algorithms etc.). 

Even there is a possibility that for the fermion determinant and for the fermion corre- 
lation functions one uses different discretizations (the first are called sea, the latter are the 
valence fermions). This is the so called mixed approach. Then the expensive, improved 
fermion is used in the valence sector, whereas for the sea fermions a faster, less improved 
is chosen. The correct continuum limit is again ensured by universality. 

The design of lattice fermion actions is hindered by the fermion doubling problem. 
Naive discretization of the continuum Dirac action yields 16 fermions on the lattice. There 
are three different ways to cure this: staggered, overlap and Wilson fermions. Let us take 
a brief look on all of them. 



Fast, but uglyQ staggered fermions 

The naive fermion determinant containing 16 fermions has an 11(4) x lf(4) symmetry, 
which can be eliminated by the staggering transformation. This reduces the degeneracy 
from 16 to four. To get one out of the remaining four fermions one applies the fourth root 
trick: 

det D [If] => detD st [lf] => (detD st [lf]) 1/4 . (1.3) 

4 Nonlocality is driven by the smallness of the quark mass. 
5 For a recent review on the staggered controversy see |5j. 
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The quark correlation functions are usuallya calculated with the four flavor operator D s t. 

The main advantage of staggered fermions is that the D s t operator has a U(l) e sym- 
metry in the massless case at any finite lattice spacing. This symmetry corresponds to the 
flavor non-singlet axial symmetries in the continuum, an important organizing principle 
in low-energy QCD. Thanks to this symmetry the spectrum of D st is bounded from below, 
making staggered algorithms well-conditioned. The bare quark mass is only multiplica- 
tively renormalized. These make the simulations fast and convenient. 

The major problem is that no local theory is known to correspond to the fourth-root 
trick, which can be a danger for the universality of the theory. Which means that it can 
happen that results of staggered lattice QCD differ from that of other discretizations. How- 
ever any attack against staggered QCD is in a hard position, it should give account for the 
remarkable agreement between staggered lattice results and real world. 

There is an explicit example ||6]|, where fourth root trick gives an incorrect result: case 
of one massless fermion. According to the anomaly the chiral condensate should acquire 
some nonvanishing value in the continuum limit with a proper fermion discretization. 
However due to the U(l) e symmetry of staggered fermions, the staggered chiral conden- 
sate is always exactly zero at any finite lattice spacing (in a finite volume), making stag- 
gered fermions fail at this setup. 

Usual calculations are done with unphysical, large pion masses, and then an extrapola- 
tion to the physical pion mass is carried out (staggered chiral perturbation theory). These 
extrapolations are controlled by several parameters (O(40) in NLO for the kaon bag con- 
stant), which can make them very ill-conditioned. 

Glorious, but slow: chiral fermions 

According to the Nielsen-Ninomya theorem eliminating the doubling problem is equiva- 
lent with violating continuum chiral symmetry ({0,75} = 0) on the lattice. The idea of 
Ginsparg and Wilson was to find a Dirac-operator satisfying 

{D, 75 } = 2D 75 D. (1.4) 

Only much later was it realized, that this relation makes possible to maintain the chiral 
symmetry at finite lattice spacing, for which an 0(a) redefinition of the chiral transforma- 
tion is needed. All continuum relations related to chiral symmetry (Ward identities, index 
theorem, low energy theorems, continuum chiral perturbation theory) are one to one ap- 
plicable at finite a using a fermion satisfying this relation. The overlap fermion and the 
fixed-point fermion are the known solutions of Eq. 11.41 whereas the domain-wall fermion 
provides an approximation to such operators. 

6 Other than staggered quark discretizations are often used in the correlation functions, these are the so 
called mixed approaches with staggered sea quarks. 



4 
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The price of these nice properties is quite high: one ends up with a non-ultralocal Dirac- 
operator (there are interactions between points at any distance), which results in 0(100) 
times or more slower algorithms compared to other discretizations. 

The determinant of a Ginsparg- Wilson Dirac-operator will have discontinuities in the 
space of gauge fields at the topological sector boundaries, just as the continuum Dirac- 
operator. It is a feature from one hand, on the other hand these jumps makes the conven- 
tional dynamical fermion algorithms with chiral Dirac-operators to slow down consider- 
ably. 

Robust: Wilson fermions 

Wilson fermions are curing the doubling problem by violating the chiral symmetry dras- 
tically. One gets several inconvenient features at a first sight: additive quark mass renor- 
malization, 0(a) lattice artefacts (the previous two discretizations have 0(a 2 )), loss of strict 
spectral bound on the Dirac-operator. Latter yields ill-conditioned, slow algorithms. The 
absence of the chiral symmetry makes necessary to evaluate renormalization constants at 
those places where it is trivial in case of staggered or overlap fermions. 

Much theoretical and numerical work was done to improve these properties: using 
smeared gauge links in the fermion action the additive renormalization can be decreased 
by two orders of magnitude, the 0(a) lattice artefacts can be reduced to 0(a 2 ) via the 
Symanzik-improvement program, the spectral bound is reported to be recovered in the 
infinite volume limit 0. 

These made possible to be competitive with the staggered discretization in speed and 
since it is in much better theoretical shape it might be the choice for the near future. 

1.2 Overview of the thesis 

The use of dynamical fermions is nowadays obligatory in lattice QCD. Developing new 
algorithms for dynamical fermions is still an active area of research. In this thesis I will 
present two dynamical fermion projects in which I have participated. 

The first topic, discussed in Chapter 2, deals with the dynamical overlap fermion project 
in detail. This chapter is partially based on the articles: 

• Z. Fodor, S.D. Katz, K.K. Szabo JHEP 0408:003,2004 

• B3 G.I. Egri, Z. Fodor, S.D. Katz, K.K. Szabo JHEP 0601:049,2006 

First I will present our dynamical overlap algorithm, then I will show how the naive algo- 
rithm fails to change topological sectors. A new algorithm is proposed and tested, which 
solves the problem. Finally the topological sector changing of the new algorithm is exam- 
ined, and attempts to improve it are proposed. I am trying to give a comprehensive review 
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of the field, which also means that only a part of the results belong to me. My contributions 
are the following: 

• Writing and developing a 5000 line C program for generating overlap fermion con- 
figurations. 

• Modifying the conventional HMC algorithm to circumvent its failure at topological 
sector boundaries. 

• Improving the stepsize dependence of this algorithm. 

• Examining the tunneling behavior of this algorithm. 

• Performing simulations and measuring the topological susceptibility in two flavor 
QCD. 

The second topic is the dynamical staggered project (Chapter 3). This is a large scale 
computation of thermodynamical properties of the quark gluon plasma. I will describe 
our choice of action and the algorithmic improvements first. Then I will show our deter- 
mination on the order of the finite temperature QCD transition in continuum limit and 
with physical quark masses. Next the transition temperature in physical units is calcu- 
lated, again in the continuum limit and with physical quark masses. These results can be 
considered as final ones modulo the uncertainty in the staggered discretization. Finally I 
will present the equation of state, but there the calculations were only done with two lat- 
tice spacings, the continuum limit is missing. The chapter is partially contained in these 
articles: 

• [HI Y. Aoki, G. Endrodi, Z. Fodor, S.D. Katz, K.K. Szabo Nature 443:675-678,2006 



• [U Y. Aoki, Z. Fodor, S.D. Katz, K.K. Szabo Phys. Lett. B643:46-54,2006 

• (H Y. Aoki, Z. Fodor, S.D. Katz, K.K. Szabo JHEP 0601:089,2006 
Again not all results belong to me, my contributions are: 

• Writing and developing a 5000 line C program for generating staggered fermion con- 
figurations. 

• Performing large scale zero and finite temperature simulations. 

• Analyzing and renormalizing the data. 
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Chapter 2 

Dynamical overlap f ermions 



Chiral symmetry is one of the most important feature of the strong interaction. Lattice reg- 
ularization and chiral symmetry were contradictious concepts for many years. Fermionic 
operators satisfying the Ginsparg- Wilson relation |13| 

{D, 75} = — D75D. (2.1) 
m 

made possible to solve the chirality problem of four-dimensional QCD at finite lattice spac- 
ing m S3 M iia - 

Several numerical studies with exact chirality operators were done in the quenched 
approximation Ifl8l [T9l |20l . The results were really compelling, but people were forced 
to work with nonlocal Dirac-operators. This made the algorithms more complicated and 
slowed them down by large factors. At the same time one could reach rather small quark 
masses, which was unimaginable with Wilson-type discretizations before. 

The life becomes even more complicated when introducing dynamical fermions with 
exact chirality. This chapter is devoted to this problem. We will going to work with overlap 
fermions ffZfl I22 ||, it is an explicit solution of the Ginsparg- Wilson relation. The other type 
of solutions, the fixed-point Dirac operators 1123 II are defined via a recursive equation, they 
are considerably harder to implement 

2.1 Overlap Dirac operator 

First we fix our notations. The massless Neuberger-Dirac operator (or overlap operator) D 
can be written as 

D = m [l + 75Sgn(H w )] / (2.2) 
1 For dynamical simulation of an approximate fixed-point Dirac operator see [24J. 
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This D operator satisfies Eq. (|2.1[) . % is the hermitian Dirac operator, % = 7sDw, which 



is built from the massive Wilson-Dirac operator, Dpy, defined by 

[D w ]xy = (4:-m )5 xy --Y^lu }l (x)(l + j }l )5 Xi y- }l + Ul(y)(l-j }l )5y /X+}l 

F 

One fermion is obtained in the continuum limit, if thq takes any value between and 2. In 
a finite volume one should be careful, that the physical branch of the spectrum of the 
operator is to be projected to the physical part of the overlap circle. 
The mass is introduced in the overlap operator by 

m 

D(m) = (!-^) D + ™- 

Sometimes it is useful to consider the hermitian version of the overlap operator. Let us 
review its properties. The massless hermitian overlap operator is H = 75D = mo (75 + 
sgn(H w )). The eigenvalues are real, the eigenvectors (|A)) are orthogonal and span the 
whole space. Due to the Ginsparg- Wilson relation 

{H, 75 } = — H 2 , (2.3) 
m 

the matrix elements of 75 satisfy: 

(A|7 5 |fi)(A + fO = — <V (2.4) 
mo 

From this equation it follows, that the zeromodes can be chosen to be chiral, furthermore 
eigenvectors with eigenvalues ±2mo have positive /negative chirality. The difference in 
the number of left and right handed zeromodes is proportional to the trace of H: 



1 ST — 

A 2^0 



n -™ = E 2^ = E^slA) " (»(+0) - n(-0)) = n(-0) - n(+0). (2.5) 



This difference is called the index of H, it can be considered as the definition of the topo- 
logical charge on the lattice. This is supported by the fact, that in the continuum limit the 
density trH converges to ~ FF. Since the right hand side can take only integer values, we 
can immediately conclude that the overlap operator cannot be continuous function of the 
gauge fields. It should be nonanalytic on the boundaries of topological sectors. 
Further eigenvalues are always coming in pairs, the vector 

'- a ^( 1 -(£f)" ,/2 (^-^)' a > < 2 - 6 > 
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is an eigenvector with eigenvalue —A. The 75 matrix leaves the subspace {|A), | — A)} 
invariant, it can be written as 



75 



/ 



A 
2m 

A 2 
(2m y 



1/2 



(2m )2 

A 
2m 



L/2\ 



/ 



(2.7) 



Since 75 is traceless on the subspace (|A), | — A)} with A 7^ {0, ±2rao}, it should be also 
traceless for the A = {0, ±2mo} space. This requires that the numbers of ±2mo eigenmodes 
satisfy 



n(+2m ) - n(-2m ) = n(-0) - n(+0). 



(2.8) 



The overlap operator squared H 2 (tn) commutes with 75. This is a trivial fact in the 
A = {0,2rao} subspaces, where even [H(m), 75] = is true. In (|A), | — A)} subspace it is 
proportional to the identity matrix 



H 2 (m) 



(2m )2 





A 2 + m 2 







(2m ) 



A 2 + m 2 



(2.9) 



2.1.1 Numerical implementation 

In the sign function of Eq. (|2.2|) one uses sgn(H w ) = H w / J H 2 ^. We usually need the 
action of the sgn(H^) operator on a given vector, which was studied many times in the 
literature [|25ll26| . The common in all algorithms is that their speed is proportional to the 
inverse condition number of the matrix H\aj. To make the algorithms better conditioned, 
one can project out the few low-lying eigenmodes of the matrix Hw and calculate the sgn 
operator in this space exactly: 

sgn(H w ) = £sgn(s)P s + sgn(QH w ), (2.10) 

where P s is the projector to the s eigenspace of H^v and Q = 1 — P s . The projections were 
done by the ARPACK code. To speed up the projections we preconditioned the problem 
with a Chebyshev-polynomial transformation [27|. We have taken the n-th order approxi- 
mation of the tanh(80(x + 1)) — 1 function in [—1, 1] interval: 

T n (x) «tanh(80(x + l)) - 1. (2.11) 

This blows up very fast around —1. If we concentrate the interesting part of the spectrum 
of Hw there, then we make the job of the eigenvector projecting algorithm considerably 
easier (its speed usually depends on the distance between consecutive eigenmodes). That 
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is we were calculating the eigenvectors of T n (H^/s^ ax — 1) instead of those of %. Since 
this is only a polynomial transformation, only the eigenvalues are different, the eigenvec- 
tors should be the same. The result is that the problem is much better conditioned using 
Chebyshev-polynomials. The speed gain was almost an order of magnitude. 

For the rest of the sign function (sgn(QHw)) one can take her/his favorite approxima- 
tion: c(QHw). We have considered two of them: Zolotarev rational function and Cheby- 
shev polynomial. 

The n th order Zolotarev optimal rational approximation for 1 / ^Jx in some interval 
[xmirv *max] can be expressed by elliptic functions (see e.g. 1128 II ), the coefficients can be 
also determined by a Remes-algorithm. A particularly useful form of the approximation 
for the sign function is given by the sum of partial fractions 



in usual cases n ~ O(10). To get the approximation of the sgn of a matrix A one has to plug 
the A into Eq. 12.121 sgn(A) ~ CzolC^)- Th e inversions appearing in this approximation all 
contain the same matrix but with different shifts (fc,). There are so-called multishift Krylov- 
space methods ||29H30f , which can solve the system of n inversions with the same number 
of matrix multiplications as needed to solve only one system. Since matrix multiplications 
dominate such algorithms, this means that practically for the cost of one inversion one ob- 
tains the solutions of n systems. When inverting the multishift system it might be desirable 
to project out even more eigenvectors than before in Eq. 12.101 and calculate the inverse in 
this subspace exactly: (s 2 + bj)~ 1 P s . 

The Chebyshev-polynomial approximation (<7"cseb) i n principle performs similarly as 
the Zolotarev rational function. Here one has to make multiplication with the Wilson ma- 
trix many times (~ 0(100)), there is no need for global summations as in Zolotarev case. 
There are architectures (eg. Graphical Processing Units I13TI ), where global summation is 
a bottleneck, it should be avoided everywhere if possible. In these cases the Chebyshev- 
approximation is a good choice. 



Hybrid Monte-Carlo (HMC, II32II ) is the most popular method for simulating dynamical 
fermions. There are many other choices possible, but HMC outperforms from all of these. 
We will base our work on the HMC algorithm. One would like to generate gauge configu- 
rations via a Monte-Carlo update with the following weight (see ll.l)) : 




(2.12) 



2.2 Hybrid Monte-Carlo 



w[U] = exp(— S 



gauge 



) det(D + D), 



(2.13) 
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where we have D + D instead of D 2 for the weight of two fermion flavors. This substitu- 
tion is legal, since the one flavor fermion determinant is real. The standard procedure to 
implement the fermion determinant is to rewrite it using bosonic fields, so-called pseud- 
ofermions ((p(x)): 

w[U] = exp(-S gauge ) J [d(p f d(p} exp(-^ + (D t D)-» (2.14) 

Now let us consider the following steps: 

1. Choose Gaussian distributed momenta P«(x). 

2. Choose a <p field according to the distribution exp(— <p^ [D^D)~ l <p). 

3. At a fixed (p background evolve the H ?( (x) gauge fields and P ? , (x) momenta using the 
equations of motion derived from the Hamiltonian 

H = I? 2 + Sg au ge + ^ + (D + D)->, (2.15) 

and from the structure of the manifold. The evolution from (U,P) fields to some 
(IZ 7 , P') via the equations of motion is usually called a trajectory. 

Iterating steps from 1. to 3. one obtains a chain of gauge configurations with the dis- 
tribution Eq. 12.131 Along a trajectory the energy and the area are conserved, moreover 
these trajectories are exactly reversible. To determine the trajectories is an infinitely hard 
problem, it requires an exact solution of the equations of motion. 

However an approximate solution of these equations can be also used to build a chain 
of gauge configurations with the exact w[U] distribution. One just has to find an area pre- 
serving and reversible integrator and integrate the equations approximately with it. With 
such an integrator at hand one trajectory will be generated via the following procedure: 

1. As before. 

2. As before. 

3. Integrate the equations of motion with an approximate integrator. 

4. Calculate the energy difference: AH = H(U' ', P') — H(U,P) and accept/reject the 
configuration with the probability min(l,exp( — AH)). 

This is one iteration of the HMC algorithm. This iteration also produces configurations 
with the same equilibrium distribution (Eq. I2.13|) as the previous one. 
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The easiest area conserving and reversible integrator is the leapfrog. The leapfrog inte- 
gration consist of making 1 / e of the following step s^: 

U(e/2)V(e)U(e/2), (2.16) 

where U operator evolves the gauge fields with the actual (fixed) momenta, wheres V 
operator evolves the momenta using the force calculated at the actual (fixed) gauge field. 
The e is called the stepsize, obviously in the e — > limit one arrives to the exact solution 
of the equations of motion. The energy is violated by 0(e 3 ) by making one leapfrog step, 
during 1/e steps of a unit length trajectory this error grows up to 0(e 2 ). The average 

energy conservation violation can be shown to be (AH) = Ce 4 + Since C is most 

generally proportional to the volume, the stepsize should be decreased as e ~ V- 1 ^ to 
keep constant acceptance. 

Using area conservation and reversibility one can show that the condition 

(exp(-AH)) = 1 (2.17) 

should be satisfied, this relation provides an easy check of the consistency of the algo- 
rithm. This condition also shows us that if one has large energy conservation violations 
(| AH | > 1), then the acceptance should be small. In this case one should decrease the step- 
size to obtain a good acceptance. The optimal acceptance rate depends on the type of the 
integrator, in case of the leapfrog and its variants the optimum is around 80%. 

Even one can drop away the area preservation property for the evolution of the gauge 
fields and momentum, in this case one has to include the Jacobian of the mapping into the 
accept/ reject step. 



2.2.1 HMC for two flavors of overlap f ermions 

In case of the overlap fermion the gauge field and momentum evolution is the following: 

U{e) : U -> exp(eP)H, V(e) : P -> P - eA 

where in the force term the A operator projects onto traceless, antihermitian matrices (in 
color indices). The complication arises in the derivative of S p f, which schematically can be 
written as 

SSpf = -xpU(D^D)xp with i/ 7 =(D + D)~V- (2.19) 

The inversion of the fermion operator tp = (D + D) -1 ^ is done by n conjugate gradient 
steps ("outer inversion"). Note, however, that each step in this procedure needs the cal- 
culation of (D^D)cp. The operator D contains cr(%), which is given for example by the 

2 In case of unit length trajectories. 

3 For studies with inverters other than conjugate gradient see l33l . 



gauge 



+ s pf ) 



(2.18) 
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partial fraction expansion (see Eq. I2.12|) . Thus, at each "outer" conjugate gradient step one 
needs the inversion of the H 2 ^ matrix ("inner inversion"). This nested type of the inver- 
sions is the price one has to pay for an exactly chiral Dirac-operator, in other formulations 
one only has one matrix inversion per force calculation. No method is known to avoid the 
nested inversions. 

The <!>(D + D(m)) = S(H(m) 2 ) derivative has a complicated form. Let us assume, that 
we treat the sgn function in the overlap operator as in Eq. 12.101 Then the derivative of 
H(m) will have two parts 

SH(m) = SiH(m) + 5 2 H{m). (2.20) 

The first part will contain the derivative of the projectors (P s = |s) (s|), this term comes not 
only as a derivative of yjsgn(s)P s , but the projectors are also there in the Qsgn(H w ) term. 
All together one obtains 

5 x H(m) = (m - m/2) £<SP s (sgn(s) - a(H w )). (2.21) 

s 

The derivative of the projector can be derived with the tools of quantum mechanical per- 
turbation theory 11341 : 

/ 1 - P, 1 - P s \ 

m = -( -ft SH W P S + P s m w - °- . (2.22) 

\H W - s H w — s J 

As we can see each projected mode brings an extra inversion of the (shifted) Wilson matrix. 
It might be safe to treat the few lowest lying modes of Hyy this way, but it is meaningless 
to calculate exactly the force of modes from the bulk of the spectrum. 

The second term in SH(m) will be the derivative of the sgn function approximation 
(52H{m)). In case of the Zolotarev approximation one has 



(2.23) 



therefore the contribution to the derivative of H(m) is 

5 2 H{m) = (m - m/2)QScr Zol (H w ). (2.24) 

In case of the polynomial approximation the formula is more complicated. 

The sgn function is needed as three different places during the HMC trajectory: in the 
force calculation one needs H(tn)~ 2 and SH(m) and in the action calculation one needs 
H(m)~ 2 again. The presence of the accept/reject step allows one to use different approxi- 
mations at different places. One can speed up the algorithm with large factors by carefully 
choosing and tuning the approximations. We were using the Chebyshev polynomial ap- 
proximation in the inversions with O(20) projected modes whereas for the SH(m) we chose 
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the Zolotarev rational approximation with 2 projected modes. The relative precision was 
always set to 10 everywhere. 

For somewhat different implementations of the standard HMC and for various im- 
provement techniques we refer here to the work of two other groups Il35ll36ll37ll38ll . 



2.2.2 Criticism of the HMC for one fermion flavor 

The HMC described above works only for a positive definite fermion matrix. This is suit- 
able for two flavors. For one flavor one can take the square root of the squared operator 
(RHMC algorithm). A different way to get the square root of the fermion determinant is 
to exploit the exact chiral symmetry of the Dirac-operator ||39l |40l . As we have seen be- 
fore H(m) 2 has n(+0) + n(— 0) chiral modes with eigenvalue m 2 , n(2mo) + n(— 2mo) chi- 
ral modes with eigenvalue (2mo) 2 and all the other modes are doublets with eigenvalue 

(•"" ~~ (2mo) 2 ) ^ 2 m2 ' ^ * S a conven ti° na l wisdom that only those configurations contribute 
to the path integral where either n(+0) = or n(— 0) = 0. If this is true then we can write 
the two flavor fermion determinant as 



n(±0) = : detH(m) 2 = det±H(m) 2 det T H(m) 2 = [det±H(m) 2 ] 2 f ^) 2 ) 



»(=F0) 



where the det± determinants have to be restricted to positive /negative chirality subspaces. 
The numerical factor at the end of the formula takes into account that the zeromodes and 
±2mo modes are not coming in chirality pairs (see Eq. |2.8|) . At this point it is easy to 
perform the square root 

/ m \»(=F0) 

detH(m) = {j^Pj) det±H(ra) 2 , (2.25) 

where the sign depends on the chirality of the zero modes of H. Since H(m) 2 is positive 
definite even on the definite chirality subspaces, there is no obstacle to introduce pseud- 
ofermions for its determinant. The contribution of the zeromodes and ±2mo modes can be 
taken into account by reweighting the observables with them. 

One has to face with the following problem. Let us consider a trajectory which starts 
with a gauge configuration where n(+0) = and ends where n(— 0) = 0. The pseud- 
ofermion is generated at the beginning of the trajectory according to the distribution: 

exp(-(p^H{m)- 2 P + (p), (2.26) 

where P+ projects on positive chirality. If one consider the reversed trajectory (starting 
from the sector, where n(— 0) =0) the pseudofermion distribution is 

exp(-(/> + H(ra)- 2 P_(/>), (2.27) 
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with P- negative chirality projector. This means that the reversed trajectory uses a different 
pseudofermion distribution than the original one. For the proof of the detailed balance 
one needs to have the same pseudofermion distribution on both ends. This algorithm 
yields the violation of the detailed balance for trajectories where the ends are in topological 
sectors with different signs. The usual way out is that the trajectories are constrained to 
that part of the phase space where e.g. n(+0) = is always satisfied. However this choice 
opens the way of a possible ergodicity breaking, which is also hard to keep under control. 



In the previous section we have shown how to set up the traditional HMC for overlap 
fermions. Performing simulations on rather small (6 4 ) lattices the acceptance rate was al- 
most zero. The strange thing was that decreasing the stepsize of the integrator did not 
help at all. Tracing down the problem, one finds that there are sudden jumps in the mi- 
crocanonical energy during the trajectories. These jumps are usually in the order of O(10) 
or larger, the trajectories where the energy violations are of this size are practically never 
accepted in the final accept/ reject step. 

These jumps occur at the discontinuity of the overlap operator, that is at the topological 
sector boundaries. The phenomena can be nicely observed in the spectrum of the hermitian 
Wilson-Dirac operator. The topological charge can be written as 



where s's are the real eigenvalues of Hpy- The charge changes when an eigenvalue of % 
crosses zero. At this point most presumably the overlap operator itself is discontinuous, 
since its trace is discontinuous. This means that the pseudofermion action also has a dis- 
continuity, which means a Dirac-delta in the fermion force. Obviously a finite-stepsize 
integrator will never notice the presence of a Dirac-delta in the force. Without any correc- 
tion one will end up with an energy violation which is roughly the discontinuity in the 
fermion action. 

One can improve on this situation. This feature is already present in a classical one- 
dimensional motion of a point-particle in a step function potential. During the integration 
one should check whether the particle moved from one side to the other one of the step 
function. If it is necessary, one corrects its momentum and position. This correction has to 
be done also in the case of the overlap fermion. The microcanonical energy, 



has a step function type non-analyticity on the the zero eigenvalue surfaces of the % 
operator in the space of link variables coming from the pseudofermion action. When the 



2.3 Reflection/refraction 




(2.28) 



H=-P 2 + S gaug e[L7] + S pi [U /( f>] = -P 2 + S[U,<p] 



(2.29) 
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When 


New momenta 


Refraction 


(N,P) 2 > 2AS 


P-N(N,P) +N{N,P)y/l -2AS/(N,P) 2 


Reflection 


(N,P)' Z < 2AS 


P-2N(N,P) 



Table 2.1: Refraction and reflection can happen to the system when approaches a zero 
eigenvalue surface of The conditions and the new momenta are indicated. P is the 
momentum before the refraction/ reflection. 

microcanonical trajectory reaches one of these surfaces, we expect either reflection or re- 
fraction. If the momentum component, orthogonal to the zero eigenvalue surface, is large 
enough to compensate the change of the action between the two sides of the singularity 
(AS) then refraction should happen, otherwise the trajectory should reflect off the singular- 
ity surface. Other components of the momenta are unaffected. The anti-hermitian normal 
vector (N) of the zero eigenvalue surface can be expressed with the help of the gauge 
derivative as 

N~(s\A(u^f)\s), (2.30) 

where A projects to the antihermitian, traceless matrices in color space. Table |2J] summa- 
rizes the conditions of refraction and reflection and the new momenta. 

2.3.1 Modified leapfrog 

We have to modify the standard leap-frog integration of the equations of motion in or- 
der to take into account reflection and refraction. This can be done in the following way. 
The standard leap-frog consists of three steps: an update of the links with stepsize e/2, 
an update of the momenta with e and finally another update of the links, using the new 
momenta, again with e/2, where e is the stepsize of the integration. The system can only 
reach the zero eigenvalue surface during the update of the links. We have to identify the 
step in which this happens. After identifying the step in which the zero eigenvalue surface 
is reached, we have to replace it with the following three steps: 

1. Update the links with e c , so that we reach exactly the zero eigenvalue surface. e c can 
be determined with the help of N. 

2. Modify the momenta according to Table |2~T1 

3. Update the links using the new momenta, with stepsize e/2 — e c . 

This means that in leapfrog step of Eq. |2.16| we have to substitute the appropriate U(e/2) 
operator with 

^mod(e/2) = U{e c )KU{e/2 - e c ), (2.31) 
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Figure 2.1: Comparison of the unmodified (triangles) and modified (boxes) leapfrogs. Up- 
per part is the energy, lower part is the lowest mode of H^y. 



where 1Z is the refraction/reflection operator, which changes the momenta according to 
Tab. 12.11 This procedure is trivially reversible and it also preserves the integration measure 
as shown in the appendix of this chapter. 

where P' = TIP is the reflected /refracted momentum. F = A ^^fr J is f ne force 
evaluated at the topological sector boundary on the starting side, which means that the 
undefined sgn(O) in F is interpreted as sgn(s— ). F' is also evaluated on the boundary, but 
with sgn(O) = sgn(s+). Obviously F' = F in case of a reflection. The energy violation of 
Eq. 12.321 can be a serious problem, since a (P,F) type quantity is in general proportional 
to the volume. Which means that one has to decrease the stepsize with e ~ V" 1 to keep 
the acceptance constant. This is much worse than the original e ~ V 1 ^ scaling of the 
leapfrog. 

Fig. 12.11 compares the evolution of the energy and lowest A for the usual and for the 
modified leapfrog. In the unmodified case there is a huge energy jump at the crossing, the 
trajectory is most probably rejected. Whereas in the modified case a reflection happens, 
and Ti. is much better conserved. There is one bottleneck however. If we correct only for 
the discontinuity in the potential as above, then the microcanonical energy violation will 
be proportional to the stepsize. That is for a U mo d(e / 2)V (e / 2) step the energy violation is 

AH = e c ((P / ,F / ) - (P,F)) + 0(e 2 ), (2.32) 



17 



2.3. Reflection /refraction 



2.3.2 Improving the modified leapfrog 

There are several ways out of the problem. We will consider two variants in this section: 

• Reflection /refraction step with Ce energy violation, the C coefficient being only an 
0(1) size number instead of 0(V), which means that the bad scaling is eliminated. 

• Reflection procedure with 0(e 2 ) energy violation. 

In the literature there exists a reflection /refraction procedure with 0(e 2 ) energy violation 
II37II , however it is somewhat more complicated than these two improvements. 

The large energy conservation violation of U mo ^ arises, since not the appropriate mo- 
mentum is reflected/ refracted. If one inserts two extra momentum updates into Eq. 12.311 

U{e c )V{e c ) • K ■ V(e/2 - e c )U{e/2 - e c ), (2.33) 

then the energy is conserved upto 0(e 2 ). This is due to the fact that both U{e)V{e) and 
T > (e)U(e) conserves the energy upto 0(e 2 ) and 1Z conserves the energy exactly. This step 
however violates the area conservation upto 0(e), there is no free lunch. 

The first idea is based on the fact, that if one makes the extra momentum updates only 
in the space, which is orthogonal to the normalvector N 

U(e c )V ± (e c ) • K ■ V ± (e/2 - e c )U{e/2 - e c ), (2.34) 

then the area is exactly conserved again (as shown in the appendix). Now the energy is 
still violated by 0(e) terms, but since the update in the orthogonal space is done 0(e 2 ) 
correctly, their coefficients is only an 0(1) number. The stepsize should be decreased only 

as e ~ V~ 1/2 . 

The second idea applies only for the reflection. It uses the following observation. In a 
one dimensional case if the time required to reach the boundary (e c ) and the time which 
is required to step away from the boundary (el 2 — e c ) are the same, then as a result of the 
reflection step, the trajectory has been exactly reversed. The area conservation, reversibility 
are obviously preserved, the energy is conserved exactly. We can try whether these remain 
true in arbitrary dimensions (the area conservation is proven in the appendix, the exact 
reversibility and energy conservation upto 0(e 2 ) are obvious). Then we should insert a 

U(e c )V(e c )lZV(e c )U(e c ) (2.35) 

step into the chain of leapfrogs, when the boundary is hit. In Eq. |2.16| we have written the 
elementary leapfrog step in the UPU form, now let us write it in the PUP order: 

V(e/2)U(e)V(e/2), (2.36) 

Now we split the evolution of the links into two parts: 

V(e/2)U(e/2) -U(e/2)V(e/2). (2.37) 
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Consider that the boundary would be crossed during one of the evolutions of the links in 
Eq. 12.371 Then replace the original leapfrog with the following: 

V(e/2)U(e/2) ■ U(e c )P(e c )HV(e c )U(e c ) ■ U(e/2)V(e/2). 

Now e c is the time to reach the boundary surface measured from the midpoint of the 
leapfrog. Thus if the crossing would happen in the first evolution then e c < 0, if in the 
second, then e c > 0. 

2.3.3 Tracing the evolution of low lying eigenmodes 

We have to trace the evolution of the low lying eigenmodes, since we are looking for the 
moment when one eigenvalue crosses zero. The eigenvectors and eigenvalues are available 
at discrete times only (once or twice per time step), therefore one has to pair the eigenvec- 
tors at time t and time t + e. We calculated the scalar products (s'(t + e)\s(t)) after each 
link update, and our recipe was the following: the s(t) has evolved to that s'(t + e) with 
which the scalar product was maximal. Of course this can break down if the time step is 
too large. It is easy to show, that to make this naive method work one has to decrease the 
stepsize as the volume is increased with ~ V . Expanding the (s'(t + e)\s(t)) one obtains: 

(s'(t + e)\s(t)) = S s , s - e(s'(t)\ d ^ 1-^ \s(t)) + 0{e 2 ), (2.38) 

at riyj — S 

where the derivative of % is simply 

Clearly the 0(e) term is proportional to the volume, therefore the e ~ V~ l relation should 
hold to be able to keep track of the evolution of the eigenvectors. If we use the derivatives 
of the eigenvectors, then we can get a considerably better scaling. That is, instead of (s'(f + 
e)\s(t)) we calculate the scalar products at time t + e/2 using the eigenvectors and their 
derivatives at t and t + e: 

(s'(t + e/2)\s(t + e/2)) = 
((s>(t + e)\-e/2±(s>(t + e)\^ (|s(f)> +e/2^|s'(f)>) +0(e 2 ). 

We have not used this formula in practice, it is expensive (to monitor N low lying eigen- 
modes one has to make N(N — 1) Wilson-matrix inversions to apply the above formula). 
Fortunately in all cases our stepsizes were always small enough that the eigenvector iden- 
tification with the naive procedure was no problem. 

The stepsize should be also small to avoid the crossing of two or more eigenvalues in a 
microcanonical time step. This happened very rarely, and since the energy violation were 
usually very large in these cases, these configurations were simply rejected. 
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2.4 Numerical simulations 



In this section we will detail the particular implementation of our overlap HMC variant. 
Afterwards we will show results on the topological susceptibility as obtained from these 
simulations. 



Gauge action 

For testing purposes the standard Wilson action was chosen as gauge action, later on we 
moved to the tree-level Symanzik-improved action. Apart from decreasing the scaling vi- 
olations in the gauge sector, the improvement is beneficial from the overlap operator point 
of view, too. Experience in the quenched case shows that improved gauge actions can dras- 
tically reduce the eigenvalue density of the negative mass hermitian Wilson-Dirac operator 
||4"Tf. Since this operator is the kernel of the overlap operator, gauge action improvement 
speeds up the overlap inversion algorithms. Since the topological sector change happens 
when an eigenvalue of the Wilson-Dirac operator changes zero, improvement reduces the 
tunneling events at the same time. Therefore one has to be careful not to overimprove the 
gauge action (as it is done for the DBW2 action). The Symanzik tree level improved action 
is the simplest improved gauge action. 



Fermion action 

The fermion action is of two overlap fermions with standard Wilson-operator as a kernel. 
After test runs with thin-links, we started smearing the links in the kernel operator via 
stout smearing procedure. The smearing reduces the fluctuations of the gauge configura- 
tion, which is again helps reducing the density of zeromodes of the Wilson-operator ||42|. 



It is an 0(a 2 ) redefinition of the gauge fields, so keeping the smearing recipe constant as 
the lattice spacing goes zero will not change the continuum limit of the theory. The stout 
smearing has the particular advantage compared to other smearing techniques, that it is 
an analytic function of the thin gauge field [43|. Therefore its derivative (which is needed 
to obtain the HMC force) can be calculated exactly. In our simulations we were using two 
levels of stout smearing with smearing parameter p = 0.15. The speedup was almost an 
order of magnitude compared to the unsmeared case. 

The negative mass of the Wilson-kernel was chosen to be — rag = —1.3 in the smeared link 
case. For smaller mo values would have been no small eigenvalues of the overlap operator, 
the topology of gauge fields would have been always trivial. The —mo = —1.3 was chosen 
to be from a small 8 4 run, where the topological susceptibility started increasing from its 
zero value at mo = and reached its plateau value around —mo ~ —1.3. 
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Algorithm 



We have tried several variants of the HMC algorithm, which were discussed in the previ- 
ous section. For all of them we were using the reflection /refraction modification in some 
way. In addition to the standard consistency tests (reversibility of the trajectories, e 2 scal- 
ing of the action and (exp(— AH)) = 1) we performed a brute force approach on 2 2 and 4 4 
lattices. We generated quenched configurations, then we explicitely calculated the deter- 
minants of H(m). These determinants were used in an additional Metropolis accept/ reject 
step. The hybrid Monte-Carlo results agree completely with those of the brute force ap- 



Now let us take a closer look on the results obtained with standard Wilson gauge action 
with standard Wilson fermion kernel in the overlap operator (results with improved action 
will be discussed later). On 4 • 6 3 lattices there is a sharp increase in the Polyakov loop 
at /3 = 5.7 (see Fig. I2.2[) , which can give a hint on the lattice spacing, since the finite 
temperature transition is usually around 200 MeV temperature (T = 1/ (4a)). This value of 
the coupling was used for measuring the topology on 6 4 lattices, which were considered as 
zero temperature lattices. The negative quark mass was set to thq = 1.6, the bare fermion 
mass was in the range m = 0.1. .1.15, the stepsize was e = 0.025 in average. Using the 
conventional HMC one would have no acceptance at all on these lattices, but modifying 
the leapfrog step according to the previous section the acceptance becomes > 70% for 
these stepsizes. At each bare mass roughly 800 trajectories were generated. The results 
are plotted on Fig. 12.31 The left panel shows the charge history. The average topological 
charge is consistent with zero for the total mass range (middle panel). 12 • 6 3 lattices were 
used to fix the scale using yq from Wilson-loops. The result is a ~ 0.25 fm for small masses. 
Pion masses were also measured and m 2 = Am with A ~ 1 is found in lattice units. Using 
the scale and the pion mass, it is possible to get the topological susceptibility in physical 
units (right panel of Fig. I2.3|) . x(m) tends to zero for small quark masses. One can compare 
these results with the continuum expectation in the chiral limit (solid line of the figure): 



proach. 



Results 




(2.40) 
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Figure 2.2: The f> dependence (right panel) of the Polyakov-loop on 4 • 6 3 lattices at m = 
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2.5 Topological sector changing 

In the previous two sections we have described a HMC algorithm for overlap fermions. 
The nonanalytic behavior of the overlap-operator at topological sector boundaries requires 
non- trivial modification of the original HMC. Our modification is able to handle the dis- 
continuity problem as shown in the numerical results section. 

As we have started simulating even larger volumes (8 4 ) with the modified algorithm, 
we had to face a new problem. In majority of the reflection /refraction steps reflection 
happened, which means that the trajectories were confined to a given topological sector 
for long times. This dramatic increase of the autocorrelation time of the topological charge 
makes the measurement of the topological susceptibility very hard and effectively also 
means the violation of the ergodicity. 

One can come up with the solution to let the trajectories continue their evolution as 
if the discontinuity in the action would not be present. This algorithm would obviously 
allow the system to tunnel between topological sectors. The price is that at the end of each 
trajectory one has to keep the exp(— AH) factors to finally reweight the configurations 
with them. The energy conservation violation is dominantly coming from the sum of the 
discontinuities in the action AH = AS a + 0(e) along a trajectory. If the system moves 
in a fixed potential, then AH will take positive and negative values equal times, since some 
times the system goes up sometimes goes down the same discontinuity. The reweighting 
would have the following form: 

r[U n ] = exp(- £ AH,), (A) = ^A[Ujr[U n ] _ ^ 

If AH/s are roughly equal times positive and negative, then the reweighting works well: 
the configurations have nearly the same weight. Unfortunately this turned out to be not 
true for the overlap HMC case, almost all AH's were positive, the configurations were 
becoming unimportant very fast in the sum of Eq. 12.411 

How could this happen? The answer is that the evolution of the trajectories is not 
done in a fixed fermion potential S ex act = — log det H(m) 2 , but in a pseudofermionic one 
Spf = (p^H(tn)~ 2 (p. The pseudofermion is not fixed, it is regenerated at the beginning of 
each trajectory. Let us take a closer look on how the pseudofermions approximate the 
fermion determinant. This will help us to understand the slowing down of the tunneling 
between topological sectors. In particular, we show that the jump in the pseudofermionic 
action overestimates AS ex act- 

Let us assume that the trajectory crosses the boundary. Let H_ and H + be the overlap 
operator evaluated on the two sides of the boundary right before and after the crossing, 
respectively. Clearly H_ and H + contain the same gauge configuration, but they differ, 
since one eigenvalue of % changes sign on the boundary. In the HMC algorithm one 
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chooses the pseudofermion field as 

<p = H-rj, (p f = rfH-, 

where tj, yf are random vectors with Gaussian distribution, in order to generate <p, cp f with 
the correct distribution. (In a real simulation one chooses new pseudofermion configura- 
tions only at the beginning of each trajectory, but for simplicity let's consider, that (p and 
<p^ are refreshed when hitting the boundary.) The jump of the pseudofermionic action now 
reads: 

AS pf = S pf+ - S pf _ = t] f (H-H+ z H- - l)rj 

The relation between AS exac t and AS p f can be obtained by the following straightforward 
calculation: 

e -AS exact = detHj = f[d V +][d V ]e-v%-rt H - H ? H --V>i = 
detH 2 j[d^)[dr])e-n^ 

= ^- 7 +(H-H+ 2 H_-l) 7 ^ + > e -{ V \H-H+ 2 H—l)r,) viv = g -( A S pf ) 

The inequality in the second line is a consequence of the concavity of the e~ x function. So 
we conclude to: 

(ASpf) > AS 

exact- 

We can examine this relation in realistic simulations, if we take into account, that there 
is a simple relation between H + and H_. Let's denote by Ao the eigenvalue of Hw which 
crosses zero at the boundary, and by |0) the eigenvector belonging to Ao- With this notation: 

H+ = H- + c|0) (0|, 

where 

tn 

c = AsgnA m (l - = — ), 
2m o 

with AsgnAo = ±2 being the jump of sgnAg on the boundary. The expectation value of the 
discontinuity in the pseudofermionic action is: 

(ASpf) = ( V \H.H- 2 H. - 1)^ = Tr(H_H; 2 H_ - 1) = 

= Tr((l -c|0)(0|H- 1 )(l — c W+ 1 |0) <0|) - 1) = -2c(0|H" 1 |0) + c 2 (0|H" 2 |0). (2.42) 

In a similar way one can get a simple formula for the exact value of the jump on the bound- 
ary: 

,-AS exact _ detHl _ 1 _ 1 _ 1 



detH 2 . det(H" 1 H_) 2 det(l - cH" 1 ^) (0|) 2 (1 - c^H^ 1 10)) 2 ' 

(2.43) 
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Figure 2.4: The jump in the exact vs. pseudofermionic action at /3 = 4.05 and m = 0.1,0.2. 
Since the average of (n,p) 2 is around w 1, topological sector changing would happen 
considerably frequently using S ex act/ than with S p f. We also indicated the probability of 
topological sector changing with the pseudofermionic action, and an estimate on the prob- 
ability using the exact action (assuming that the two algorithms would behave the same 
way except for the boundaries). 



Eq. (|2.42|) and Eq. (|2.43|) offers a numerically fast way to determine both action jumps, 
since one needs only one inversion of the overlap operator to obtain both of them. 

For illustration we made a scatter plot (Fig. 12 .4|) from a 6 4 lattice at two different masses. 
From the joint distribution of AS eX act/ AS p f we can understand why are the tunneling events 
are so rare. Topological sector changing occurs when the HMC momentum of the system 
in direction of the topological sector boundary surface is large enough to "climb" the dis- 
continuity (see Tab. I2.1|) . The momentum squared is usually an 0(1) number. As we can 
see on Fig. 12.41 the AS p f distribution overestimates the real discontinuity AS eX act with or- 
ders of magnitude. Therefore a crossing which would be possible with AS eX act becomes 
impossible with AS p f. The HMC which uses AS p f stucks into a given topological sector. 
The overestimation becomes worse with lowering the quark mass. 

One way to cure this is to use several pseudofermion estimators instead of one ||36|. 
More pseudofermions mean smaller spread of the pseudofermionic action distribution, 
therefore the overestimation is smaller, too. However the computational time also in- 
creases with the number of extra fields. Obviously the best would be to use the exact 
action in the simulations, but only its discontinuity on the boundary can be calculated eas- 
ily (the calculation of the exact fermion determinant is an 0(V 3 ) operation in general). In 
the following two subsections we show two different ways to use the exact action jump 
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instead of its pseudofermion estimator in the simulations. Both of them are inexact, the 
errors present in the measured quantities are of 0(e 2 ). 

2.5.1 Using ASexact to sew together simulations with fixed topology 

Let us write the partition function in the form (assuming a vanishing 6 parameter): 

CO 

z= E Z Q> 

Q=-oo 

where Zq is the partition function of the topological sector Q. The expectation value of an 
observable: 

(n , _ L Q z Q (Q) Q _ got (Q)q 

where the restricted expectation value (0)q is 

(0) Q = J- / [dlZ] Q 0[U] detH Q exp(-S g ). 

For reasons which will be clear later the integration goes not only over the configurations 
with Q charge, but also over the boundary of the topological sector as well (though the 
boundary has only zero measure in this case). When calculating the partition function in 
a given topological sector the following boundary prescription is used: we define the de- 
terminant on the boundary as the limit of determinants approaching the wall from the Q 
side (det Hq). If the measurement of the quantities Zq + i / Zq would be possible, then we 
could recover Zq / Zq for any Q. With these in hand, we would need only the restricted ex- 
pectation values (0)q, whose measurement doesn't require topological sector changings. 

Measuring Z Q+1 /Z Q using AS exa ct 

Now we will show a way to measure Zq + i / Zq. It will make use of the fact, that we can 
calculate easily AS eX act on the boundary of topological sectors (see Eq. (|2.43|) ). The pseud- 
ofermionic action is only used to generate configurations in fixed topological sectors, so its 
bad distribution for the jump of the action will not effect us. (In the following formulae AS 
will automatically mean AS eX act-) The main idea is the following: an observable measured 
in sector Q is inversely proportional to Zq and an observable in Q + 1 is to Zq + \. If the 
observables in the two sectors are concentrated only to the common wall separating the 
two sectors, then from the ratio of the two expectation values one can recover the ratio of 
the two sectors. 
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First let us measure in the Q sector an operator, which is concentrated to the boundary: 

<<Wi f >Q = -f- I [dU} Q S QiQ+l F[U] detH 2 Q exp(-S g ), (2.44) 



where we introduced the distribution 5q : q + \, a Dirac-c5, which is equal to zero everywhere 
but on the Q,Q + 1 boundary. Then let us measure another operator G on the same wall 
(thus on the boundary separating sectors Q and Q + l), but now from the Q + 1 sector: 

(S QrQ+1 G) Q+1 = / [dU] Q+1 S QrQ+1 G[U] detH 2 Q+1 exp(-S g ). (2.45) 

Z Q+1 J 

The wall is the same (i.e. [dU] q8q,q+i = [dU] q+i^q+i) in both cases, however due to our 
boundary prescription the determinants are different on it. Therefore if F and G satisfies 

F[U] detH^lT] = G[U] detH^ +1 [(J] (2.46) 

for configurations on the boundary, then the ratio of Eq. (|2.44]> and Eq. (|2.45|) gives us 

(%Q+1 F )q _ Z Q+l 



(^Q,Q+1 G )<2+1 Z 



Q 



(2.47) 



Choosing F[U] and G[LT] functions 



The easiest choice is G(U) = 1 and F(U) = detHg +1 / det FZq = exp(— AS), the ratio of 
sectors becomes: 

\ d Q,Q+VQ+l 

This choice is still not optimal, since the measurement of the numerator is problematic, if 
the distribution of AS extends to negative values. The exponential function amplifies the 
small fluctuations in the negative AS region, which can destroy the whole measurement: 
a very small fraction of the configurations will dominate the result. As a consequence 
one ends up with relatively large statistical uncertainties. With a slightly different choice 
of F and G we can improve on the situation. With F(U) = ©(AS — x)exp(— AS) and 
G(U) = ©(AS — x) we can omit the problematic part of the AS distribution (the values 
smaller than x) from the measurement, and we get: 

(% Q+1 exp(-AS))£ s> * 
\ d Q,Q+VQ+l 

The price of this choice of F, G is that we do not make use of the AS < x part of our data 
set. The value of x can be tuned to minimize the statistical error. 
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Let us note that Eq. (|2.46|) can be viewed as a detailed balance condition on a given 
U configuration between Q and Q + 1 sector (F and G are just the "transition probabil- 
ities"). This can give us a hint, that the Metropolis-step is a good a solution for F, G: 
F = min(l,exp( — AS)) and G = min(l,exp(AS)). The ratio of sectors is simply: 



The inconvenient part of the distribution (AS < 0) is cut off, however in contrast to Eq. 
I2.49l all configurations are used to get the expectation values. 

Expectation value of a Dirac-delta type operator 

Let us discuss briefly that in the framework of HMC, how to measure an expectation value, 
which contains a Dirac-delta on the surface. The important observation is that one can use 
the pseudofermionic action in the HMC to get the fixed topology expectation values in Eq. 
l2.48H2.49l 12.501 Inside a topological sector the behavior of the pseudofermionic estimator is 
not an issue, we can use it instead of S ex act as usual. In practice it is not possible to measure 
an operator containing a Dirac-delta on the boundary surface on configurations generated 
by the pseudofermionic HMC, because none of them will be exactly located on it. If we 
would be able to exactly integrate the equations of motion, then all inner points of the 
trajectories could have been taken into the ensemble. Those ones also, which are located 
exactly on the surface. Here one would pick up a contribution from the Dirac-delta to the 
above expectation values, at the inner points the contribution would be zero. In the real 
case the trajectories differ by 0(e 2 ) from the exact ones^. Here using the above procedure 
(measuring the F[U] and G[U] operators on the boundary and summing them up along 
the trajectories) one makes 0(e 2 ) errors in expectation values. 

Summarizing the new technique 

We have achieved our main goal: without making expensive topological sector changes we 
can obtain the ratio of sectors (see Eq. l2.48H2.49ll2.50|) . The key point is to make simulations 
constrained to fixed topological charge, and match the results on the common boundaries 
of the sectors. Since no sector changing is required, the inconvenient distribution of the 
pseudofermionic action jump on the boundary will not effect the measurement of the ratios 
of sectors. The exact action is needed only on the boundary: the formulas 12.481 12.491 12.501 



Obviously an important issue for this new method is whether topological sectors de- 
fined by the overlap charge are path-connected or not. In [44| it has been proven that 

4 In order to have only an 0(e 2 ) difference one has to use an improved modified reflection step as de- 
scribed in the previous section. 



Z Q+l/ Z Q 




(2.50) 



(%Q+l min ( 1 ' ex P( AS )))Q+l' 



require AS. 
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Abelian lattice gauge fields satisfying the admissibility condition can be classified into con- 
nected topological sectors. No result is known for non- Abelian groups or non-admissible 
gauge fields. (Though there are some concerns on the structure of the space of non- Abelian 
lattice gauge fields ||45l .) If configurations with the same Q would not be continuously 
connectable in sector Q, then our assumption that we make measurements on the com- 
mon boundary of sectors could be violated. It could happen, that the wall sampled from 
sector Q does not coincide with the wall sampled from Q + 1. Moreover the fixed sector 
simulations would also violate ergodicity in this case. Let us note here that the large au- 
tocorrelation time for the topological charge in the conventional pseudofermionic HMC 
effectively also causes the breakdown of ergodicity. In case of non-connected sectors one 
can cure these problems by releasing the system from a sector after a certain amount of 
time and closing it to another. 

2.5.2 Using AS exa ct in R-algorithm 

In the following we will describe another technique, which uses the AS eX act and can cir- 
cumvent the critical slowing down of the topological sector change. If one does not insist 
on an exact algorithm, then an R-algorithm [|46| where the AS ex act's are taken into account 
can be a particularly good choice. Let us describe it shortly. Instead of evolving the trajec- 
tory in a pseudofermion potential (see Eq. |2.18|) , one can try to estimate the exact force by 
a random vector: 



Usually one estimator (R) per integrator step is used, so the approximation might be poor. 
If the stepsize goes to zero, then on a fixed time interval the number of estimators will 
diverge making the approximation exact. Since there is no recipe, how to make the R- 
algorithm at finite stepsize exact (like the accept/ reject step in the HMC algorithm) the 
stepsize extrapolation is a necessary ingredient. The stepsize error scales with 0(e 2 ). 
When a trajectory hits the topological boundary surface, then one just has to modify the 
trajectory according to the reflection/ refraction rules, but now one can use the AS exact dis- 
continuity instead of a badly behaving estimator (eg. AS p f). The modified leapfrog step 
is not necessarily to be an exactly area conserving one (since stepsize errors are already 
present). But still it is required, that the errors caused either in the energy or in the area 
conservation are minimal (a good candidate is the leapfrog-in leapfrog-out, which con- 
serves the energy upto 0(e 3 ) and the area upto 0(e 2 )). 




(2.51) 
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Figure 2.5: Left panel: a typical optimization procedure of the lower limit (x) on AS in 
the formula (|2.49|) . The statistical error of the ratio Z\/Zq shows a minimum as the func- 
tion of x, which is considered as the optimal value. Right panel: Bare mass dependence 
of topological susceptibility using three different methods on 6 4 lattices. The points corre- 
sponding to the same mass were slightly shifted vertically for clarity. Result based on our 
new technique and Eq. 12.491 is on the left, based on Eq. 12.501 is in the middle, the standard 
pseudofermionic HMC is on the right. The simulation parameters are from [|47|. 



2.6 Numerical simulations 2. 

In the previous section we described two methods, to solve the topological sector changing 
problem of pseudofermionic HMC simulation. We were extensively using the first one (see 
subsection I2.5.1|) . Here we describe the details of these simulations, and finally give the 
topological susceptibility in physical units measured on 8 4 and 8 3 x 16 lattices. 

Simulations were done using unit length trajectories, separated by momentum and 
pseudofermion refreshments. The system was confined to a fixed topological sector in 
each run, we reflected the trajectories whenever they reached a sector boundary. The end 
points of the trajectories obviously follow the exact distribution in a given sector, usual 
quantities can be measured on them. We compared a few observables (plaquette, size 
of the potential wall) in a given topological sector, but in different runs. We have not 
found any sign indicating that the sectors were disconnected. When calculating the ratio 
of sectors using Eq. 12.481 or Eq. 12.491 or Eq. 12.501 we integrated along the trajectories, this 
quantity will be burdened by a step size error. We carried out simulations at one stepsize. 

In case of large enough statistics the value of Zq + \ / Zq should be the same, indepen- 
dently which of the three formula 12.481 12.491 and 12.501 was used to calculate it. We omit 
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Eq. 12.481 in the following, since it is hard to give a reliable error estimate on the expecta- 
tion value of exp(— AS), if AS can be arbitrary negative number. Eq. 12.491 still measures 
exp(— AS), but with a lower limit (x) on AS. Smaller limit yields a smaller and more re- 
liable error, however the statistics is decreased at the same time. One can tune the value 
of x, so that the statistical error takes its minimum. A result of a typical optimum search 
can be seen on the left panel of Fig. 12.51 The optimal value can be compared to the one 
obtained from Eq. 12.501 On the right panel of Fig. 12.51 the two new topological suscepti- 
bilities and the one calculated by using traditional pseudofermionic HMC [47| are shown. 
The agreement is perfect. Comparing these results with those of the HMC, we conclude 
that the stepsize effect is negligible (at least at our present statistics). Let us compare the 
amount of CPU time of the two different methods for roughly the same statistical errors 
(see Fig. I2.5[) : the conventional HMC consisted 500-1000 trajectories (500 for the smallest, 
1000 for the largest mass), whereas we generated less than 200 at each mass for the new 
method. Moreover it is important to emphasize in this context that the new method can be 
efficiently parallelized. 
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Table 2.2: Topological susceptibility measured on 8 4 lattices in the second and third col- 
umn. The further columns contain the Sommer-scale, pion mass, pion mass times box size 
and number of trajectories on 8 3 x 16 lattices. 

To measure the topological susceptibility on 8 4 lattices we generated configurations 
with tree-level Symanzik improved gauge action (/3 = 4.15 gauge coupling) and 2 step 
stout smeared overlap kernel (p = 0.15 smearing parameter, the kernel was the standard 
Wilson matrix with m,Q = 1.3). We performed runs in sectors Q = ... 3 (based on the 
measured Z3 / Z2 we can conclude, that the contribution of Q > 4 sectors are small com- 
pared to statistical uncertainties). For the negatively charged sectors we used the Q — > — Q 
symmetry of the partition function. The bare masses were m = 0.03, 0.1, 0.2 and 0.3, at each 
mass approximately 1000 trajectories were collected. The average number of the topolog- 
ical sector boundary hits was around 1.5 per trajectory. We calculated the ratio of sectors 
using Eq. 12.491 and Eq. 12.501 The result for the topological susceptibility can be seen on 
Fig. 12.61 (see also Table |Z2|). It is nicely suppressed for the smallest mass. To convert it into 
physical units, we made simulations on 8 3 x 16 lattices. We measured the static potential 
by fitting the large time behavior of on and off-axis Wilson-loops. Then fitting it at interme- 
diate distances we extracted the value of Sommer-parameter. We also measured the pion 
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m (m,r )« 

Figure 2.6: Topological susceptibility as the function of quark mass on 8 4 lattices in lattice 
units (left), and in physical units as the function of pion mass (right). Scale fixing and mass 
measurements were done on 8 3 x 16. The error bars on the right plot do not contain the 
errors of scale fixing. The line is the leading order chiral behavior in the continuum. 

mass (see Table |Z2|). Since our statistics was quite small on these asymmetric lattices, the 
errors are large. Note, that in order to get the mass-dimension 4 topological susceptibility 
in physical units, one has to make very precise scale measurements. 

When interpreting the results, one should keep in mind, that the volume is small, and 
the lattice spacing is large. Note however, that smeared kernel overlap actions show nice 
scaling behavior and good locality properties already at moderate lattice spacings [42ll48l . 

2.7 Discussion 

In this chapter we have given a summary of the work to implement a dynamical overlap 
fermion algorithm. The lattice index theorem of the overlap Dirac-operator is a very nice 
feature, however it has its bottleneck. The operator is nonanalytic at the topological sector 
boundaries, which makes the conventional dynamical fermion algorithm (HMC) break 
down. We have proposed, implemented and tested a modification which is able to handle 
this nonanaliticity. Examining the properties of the modified algorithm carefully, we have 
made a few improvements on it. One of them was an improvement of the acceptance ratio, 
the other is connected to the slow topological sector changing of the algorithm. 

Even with these improvements the simulation with dynamical overlap fermions is in 
an exploratory phase. Other fermion formulations are considerably faster than the overlap. 
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There are two major problems at the moment. 

1. The first bottleneck is that the construction of the overlap operator is a very expensive 
procedure, it scales with V 2 (as one can see in Ref. Il49i but the extra V factor can 
be expected, since the number of zeromodes increases with the volume). Therefore it 
is very hard to imagine a dynamical fermion algorithm with better scaling behavior. 
As it was mentioned in the introduction the algorithms for conventional fermion 
formulations scale usually with V 5 ^. 

2. The second bottleneck is handling the nonanaliticity of the overlap operator. The 
most simple modification of the conventional HMC (as described in the chapter) can 
easily bring extra V factors in the scaling. There exists modifications improving the 
situation as we have seen, but they are really cumbersome. The problem is that the 
more sophisticated an improvement is, there are more ways to go wrong. The nice 
feature of the HMC, the robustness will be lost. 

Without a solution of the first issue at hand (which would mean to get rid of the nested 
inversion), one can simply accept that the overlap dynamical fermion algorithm will scale 
at least with V 2 . At this point new algorithms might come into play where presumably dif- 
ferent problems have to be solved. If the only gain is that one can forget the discontinuities 
in the overlap operator (second issue), it might worth to change. 

2.8 Appendix: area conservation proof 

The leapfrog is trivially an area conserving mapping in the phase space, since the increase 
of the momentum depends only on the actual coordinates, and the change in the coor- 
dinates depends only on the momentum. In case of the modified leapfrog the difficulty 
arises since e.g. in the first step the updates of the link variables are depending on the 
actual links through e c . Similarly the momentum update also depends on the momentum 
through the normal vector. 

In order to keep the discussion brief, first let us start with a Hamiltonian system in the 
N-dimensional Euclidean coordinate space. This shows the basic idea of the proof in a 
transparent way. 

We solve the equations of motion with a finite stepsize integration of the following 
Hamiltonian: 

H = -VaVa + S (sgnA%)) , 

where q a , p a (a = 1 . . . N) are the coordinates and the momenta. M depends only on the co- 
ordinates and the action S is a smooth function (note that q a , M and S are analogous to the 
links, the fermion matrix and the fermionic action, respectively). The standard leap-frog 
algorithm can be effectively applied to this system, as long as the trajectories do not cross 
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the zero eigenvalue surface of M (A(^) = 0, where A(^) is the eigenvalue with smallest 
magnituddj). 

We have to modify the leap-frog algorithm, when the coordinates reach the zero eigen- 
value surface. Instead of the original leap-frog update of the coordinates, where the con- 
stant p a momenta are used for the time e/2, we first update the coordinates with p a until 
the surface, then we change the momentum to p' a , which is used to evolve q a for the re- 
maining time. In case of refraction one has the following phase space transformation: 

q' = q + e c p + (e/2 - e c )p' (2.52) 
p' = p — n(np) + n(np') + h, 

where n is the normalvector of the surface, AS is the potential jump along the surface, and 
(np') 2 = (np) 2 — 2AS. e c is the time required to reach the surface with the incoming mo- 
menta p. h is a vector orthogonal to n and depending on q, p only through e c or quantities 
which measured on the eigenvalue surface. The h might be needed to improve the energy 
conservation of the leapfrog (see Sec. |2.3|) , eg. one can use 

h = -e c QF- - (e/2 - e c )QF+, (2.53) 

where the F± forces are measured on the eigenvalue surface with setting sgn(A(e c )) = 
sgn(A(e c ± 0)). Q a j, = b a \, — n a n\, is simply the orthogonal projector to the surface. 

First let us concentrate on the q, p dependence of e c . e c (q,p) is determined from the 
condition \(q + e c (q,p)p) = 0. One obtains the partial derivatives of e c with respect to 
q, p by expanding this zero eigenvalue condition to first order in Sq or dp. First take the dq 
variation: 



X(q a + e c p a + Sq a + ^5q b p a ) = A(q + e c p) + 

dq h dq a 



de 

(Kb + ^Pa^qb = (2.54) 



q+e c p 



dq 



b 



Since the normalvector is just 



3A 

n a = t— 

dqa 



q+e c p 



9A 
dq 



we have for the partial derivative of e c with respect to q: 

de c _ n a 
dq a (np) ' 

Similarly one gets for the partial derivative with respect to p: 

de c n a 
-e c - 



dp a (np) ' 



5 We do not deal with the possibility of degenerate zero eigenvalues which appears only on a zero measure 
subset of the zero eigenvalue surface. 
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There is an important identity between the q and p derivatives of a function, which 
depends only on q + e c (q,p)p. (Two examples are n and AS.) Let us evaluate p and q 
derivatives of an arbitrary g(q + e c {q, p)p) function: 



dq a dq b 

dg _ dg 
dp a dq b 



fx , 9e c \ d 8 

q +e c p Ma Mb 

fx d£ c \ d g 

q+e lP d Va Mb 



q+e c j> 



q+e c p 



_ n a p b , 
~ (np) h 

{up) )€c ' 



(2.55) 
(2.56) 



which gives 



M_ 



dq„ 



(2.57) 



Now we can consider the four different partial derivatives required for the Jacobian: 

/ = 



dq' dq' 

dq dp 

dp 1 dp 1 

dq dp 



whose determinant gives the change in the Euclidean measure d N qd N p due to the given 
phase space transformation. Introducing 

7 = d A 
dq b 

one incorporates all terms which arise from the q dependence of the normalvector and AS. 
In case of a straight wall with constant potential jump and h = this matrix vanishes. 
(Clearly, for QCD with overlap fermions this object is very hard to calculate; they usually 
require the diagonalization of the whole % matrix ). Using Eq. 12.571 one can recognize 
the Z matrix in the other three components of /. Denoting 



X, 



ab 



Qab + 1 



2AS 

(np) 2 



1/2 



n a n b + 



h a n b 
(np) 



y a b = Qab + 1 



2AS 

(np) 2 



-1/2 



n a n b 



(2.58) 
(2.59) 



The useful property of X and Y that the determinant of their product is 



det(XY) = det 



hb + 1 



2AS 

(np) 2 



-1/2 



h a n b 
(np) 



1+1 



2AS 

(np) 2 



-1/2 



(hn) 
(np) 



(2.60) 
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which means that it is trivial for the (/in) = case. In terms of the X, Y and Z matrices 
the Jacobian is very simple. We can split it into 2 parts: the first term contains all X and Y 
factors and has determinant one and all Z factors are in the second term: 

X e c X+{e/2-e c )Y\ , f(e/2-e c )Z (e/2-e c )e c Z 



— ~^ + ^ v _ — £ — j. (2.61) 

Let us introduce /' as the product of / and the inverse of its first term. Simple algebra gives: 

J'=(l J) ® 1 + E ® ^Y^Z, (2.62) 

where E is defined as 

\l/e c 1 

E has an eigenvector V\ oc (e C/ —1) with zero eigenvalue. The Vi oc (1, e c ) vector is orthog- 
onal to V\ and has the property to give zero in the product V2EV2 = 0. In the orthonormal 
basis given by v\ and V2 }' has the form: 

''=(£ ?)® 1+ (o ^f^^r-'Z, (2.63) 

thus det /' = 1. Since / and /' differs only in a matrix with determinant one, we arrive 

det J = l, 

thus the transformation Eq. [Z52] preserves the integration measure. 
The transformation for reflection is given by 

q' = q + e c p + (e/2 - e c )p' (2.64) 
p' = p — 2n(np) + h. 

h can be chosen as in Eq. 12.531 but now we have F_ = F+, since at reflection the sgn 
function does not change sign. One can obtain the Jacobian of reflection by simply making 
the 



2AS ^ 1/2 



[n P y 



(2.65) 



substitution in the Jacobian of the refraction (Eq. |2.61|) . Then it is easy to see that the 
det / = 1 holds for the reflection case, too. 

Finally let us consider a modified reflection, which makes only 0(e 2 ) error in the energy 
conservation (see Sec. |2.3[) . The phase space transformation can be written as: 



q' = q + e c p + e c p' (2.66) 
p' = p — 2n{np) + h, 
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The h which is needed to ensure energy conservation upto 0(e 2 ) is the following 

h = -2e c QF_. (2.67) 

This comes from Eq. 12.531 and using that the inward and outward updates now take the 
same time (e c ). h automatically satisfies (hn) = 0. The Jacobian is very similar to the 
Jacobian of the reflection procedure above (ie. the one obtained from Eq. 12.611 with the 
( ) 1/2 -> -1 substitution): 

Instead of e/2 — e c we have e c everywhere and the X matrix is substituted by Xi: 

_ (2Qp + h) a ti b 

[ X l\ab - d ab • (2-69) 

Xi has a trivial determinant detXi = 1, since (nh) = (n,Qp) = 0. From here the proof 
goes in the same way as above. One concludes to det / = — 1, where the minus signal comes 
from det Y = — 1. 

The proofs for the SU(3) cases were carried out in a completely analogous way. The 
only difference was the appearance of factors associated with the group structure of SIT (3) 
which all canceled in the final result. Thus, we conclude that the suggested modifications 
of the leap-frog conserve the integration measure. 



2.9 Appendix: Classical motion on an SU(3) manifold 

In this appendix we briefly discuss the Hamiltonian formulation of a system, which co- 
ordinates are elements of a G = SIT (3) group. In particular we will provide formulas to 
calculate the Jacobian of some map in the phase space. Some parts of the appendix closely 
follow Ref. El. 



2.9.1 Differential geometry on a Lie-group 

If the coordinates of a system are elements of a Lie-group manifold (g a G G), then T g G 
is the space of tangent vectors at point g, this is the vector space of velocities (with local 
coordinates g a ). 

Let us consider a few relevant mappings which arise due to the Lie-group structure of 
G. There is a natural mapping called the right translation 

R g :G^G h^hg, (2.70) 

6 In the previous reflection recipe, det X = — 1 was also true, so all together one ended up with det / = 1 . 
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the corresponding derivative mapping Rg*{h) : T/ 7 G — > TWG is a linear transformation 
which has the following matrix in local coordinates: 

The pullback of Rg is in certain sense going in backward direction as in the case of the 
derivative mapping, since 

R* g (h) : n g G - T h G cc(hg) - (5(h) : f}(h)(v) = oc{hg) (R g *(h)v) , (2.71) 

for all v vectors in the tangent space T/,G. Here oc and f> are 1 -forms, linear functionals 
acting on vectors. The group element dependence is indicated in the (...) parentheses, 
whereas the vector, which they act on, is in the (...) bracket. 

A vector field is right invariant, if v(hg) = Rg*(h)v(h) is fulfilled. There is a one 
to one correspondence between right invariant vector fields and the elements of the Lie- 
algebra of the group (v(g) <-> v(l) £ LA(G)), thus they are elements of a linear space. The 
Lie-bracket of two vector fields (v and w) measures the noncommutativity of two flows 
(one parameter G — > G maps, whose derivatives are the vector fields themselves). It is 
again a vector field: [v,w] = u, or in local coordinates it is u a = w\^ b v a — v b d b w a . For right 
invariant vector fields the bracket is also right invariant, thus if r& is a basis in the linear 
space of right invariant vector fields, then 

VaM = c C AB r c . (2.72) 

A 1-f orm field is right invariant, if R*ct = ol, that is oc(hg) {R g * (h)v) = oc(h) (v) for all v 
vectors in tangent space T^G. There is a one to one correspondence between right invariant 
1-f orm fields and 1 -forms over the tangent space at the unit element (cc(g) <-> ol(1), so 
that oc(l)(v) = a(g){Rg*v)). In order to prove an important identity for right invariant 
1-forms, we need a little preparation. If a = oc a (g)dg a is a 1-f orm field, then its derivative 
is da = dfrttadgi, A dg a . Its pullback corresponding to a mapping R is R*a = eif,(k)d a ki,dg a 
with the k = R(g) notation. Then 

d{R*cc) = {d d a b {k)d c k d d a k h + a b (k)d c d a k b )dg c A dg a = 

= d d a h {k)R* {dg d A dg h ) = R*doc, 

where we have used the antisymmetric property of the wedge product. Using the above 
equation it is easy to see that the derivative of a right invariant 1-form is also right invari- 
ant: 

R*da = dR*a = da. 
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This means that if we take Qa as a basis in the space of right invariant 1-forms 0, then dq^ 
should be expressible in terms of qb A qc- So let us calculate the 2-form cIqa ° n two basis 
vectors in local coordinates: 

d Q A (r B ,r c ) = d bQ A (dg b A dg a ){r B ,r c ) = (d bQ A - d aQ A )r B r c a = 
rtMQara) ~ r%{qtrl) - r B Q A d b r c a + r c aQ A d a r B . 

In parentheses we have Sac an d Sab due to the normalization, therefore only the last two 
term remains. These two gives —Q A ([r B ,r c ]), which yields the following result (Maurer- 
Cartan structure equation): 

dQ A = -\c A c Q B AQ C . (2.73) 



2.9.2 Hamiltonian dynamics 

The Lagrangian of the system is a real valued function on the tangent bundle (L : TG — > 
1Z). The derivative of the Lagrangian in the direction of the velocities is a differential form, 
which maps from TgTgG ~ TgG to the real numbers (ie. it is an element of the cotangent 
bundle T*G). Its local coordinates are j^s which are identified as the canonical momenta 
(p a ). Since the momenta are coordinates of linear forms on TG, the Hamiltonian phase 
space is the manifold T*G with local coordinates {g a , p a }- 

The T*G manifold is symplectic, ie. we have a 2-form w on T*G which has vanishing 
derivative: 

co = d\TQ a pa\ ^dcv = 0. (2.74) 

According to the Maurer-Cartan equation, the 

\ 

v = Y^Qa^dp a + -p a c a bc Q b A Q c (2.75) 
a 

relation holds. From the symplectic structure follows, that there is a one to one correspon- 
dence between vector fields (v) and 1-form fields (a): 

v «-> ol a{w) — co(v,w). 

for all w vectors. 

7 It is normalized so, that Qa(tb) = Sab is satisfied at the identity. Due to right invariance, the normaliza- 
tion will hold on the whole group. 
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The equations of motion arise through a Hamiltonian function (H) and the symplectic 
structure. The change in the Hamiltonian is described by the derivative 1-form dH. Along 
the vector field h, which corresponds to the 1-form dH through the symplectic structure, 
the Hamiltonian is conserved: 

dH(h) = co{h,h) = 0. 

In order to determine h we use the right invariant vector and 1-form basis on the group. In 
this basis the Hamiltonian vector field h and an arbitrary vector field v has the following 
form: 

d d 
h = h a r a + h a —, v = v a r a + . (2.76) 

The derivative 1-form of the Hamiltonian dH can be written as 

dH 

dH = dH(r a )Q n + —dp a , 
dp a 

where dH{r a ) is just the r a directional derivative of H. Now it is easy to see that 

dH 

co{h,v) = h a v a - h a v a + c a bc p a h h v c and dH(v) = v a dH(r a ) + v a ^ — 

oPa 

holds. Equating coefficients of v a and v a we get the result for h: 

r)H r)H r) 

h =w/" +(4 ' Vc w b - dH{r " )] w.- (277) 

The integral curve corresponding to the vector field h describes the motion of the system 
in the phase space as the time (t) goes on. The equations of motion are the differential 
equations for {g a , p a } coordinates which is solved by the integral curve: 

dH dH 

p a (t) = dp a {h) = -dH{r a ) + c c ba p c — and g a (t) = dg a {h) = ^—dg a {r b ). (2.78) 

op b oPb 

2.9.3 Volume element and phase space maps 

The integral curved corresponding to the Hamiltonian vector field preserves the symplectic 
structure g(t)* co = to, which means cv(g(t))(g(t)*v,g(t)*w) = cv(g(0))(v,w). Moreover 
higher "wedge" powers are also preserved, notably the largest one 

Q = co d = co A co A • • • A co = Qi A (> 2 • • • A Q d A dp 1 A dp 2 • • • A dp d/ (2.79) 



3 For simplicity the notation of the integral curve is g(t) instead of {g(t), p(f)}- 
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with d being the dimension of the group. O is the volume element of the phase space, it 
is the wedge product of the Haar-measure of the group and a Euclidean volume element 
(d d p). The g(t)*Ci = Q property is usually called area conservation. 

Let us consider a phase space map / : T*G — > T*G with {g, p} coordinates mapped to 
{G,P}. Since O is the only 2rf-form on T*G, f*Cl is proportional to O. The proportionality 
constant describes the change in an infinitesimal phase space volume under the map /. By 
definition 

(gp) . . . , v {2d) ) = O(GP) {f t v {1) , . . .,f*v m ), (2.80) 

with/* being the derivative mapping of/, vr^'s are arbitrary vectors in T g p(T*G) tangent 
space. In the usual basis (see Eq. I2.76|) an f*v vector can be written as 

r d 

f*v = v a f*r a + v a f* 



Eq. 12.801 is actually a 2d dimensional determinant, in which we have to deal with the 
following types of objects: 

Q a {f*v) = v h Q a {f*r h ) + v h q a {f*^—) and dP a (f*v) = v b dP a (f*r b ) +v b dP a (f* ' 



opb dpb 
Based on these relations the determinant of Eq. |2.80| is 

n(GP)(f*v {1) ,...,f*v i2d) ) =det(}v {1)l }v i 2 ) ,...,}v {2d) ) = {det})n(gp){v {1)/ . . .,v {2d) ), 
where / hypermatrix was introduced as: 

(Qa(f*r b ) <?«(/; 4>\ 

det / is the proportionality constant that we were looking for. 
2.9.4 Formulas in matrix representation 

In practice the dynamics is treated in terms of matrices instead of independent real param- 
eters. The group variables are represented by unitary matrices. A possible parametriza- 
tion is U(g) = exp(g a T a ) with T a traceless, antihermitian matrix basis. The momentum 
becomes a traceless, antihermitian matrix TL(p) = p a T a - Let us consider the r A directional 
derivative of U(g): 



dU(r A ) = ^-r A (g) = ^dg b (R g *(l)r A (l)) 



dU d(hg) b 



dg b dh c 



r A m _ dU(hg) 
h=0 an c 



r A (l) = T c Ur A (l) = T A U. (2.82) 



h=0 
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We have used the right invariance, the form of Rg* in local coordinates, the explicit form of 
U(g) and finally we have fixed the local coordinates of the r A basis at the identity (r A (1) = 
5 Ac)- We will also need an equation similar to the above 



dg^_ 



d(hg) b dh A 



dU?(hg) d(hg) h 
dUQi) 



dh A 



h=0 

dh A 



dU?(hg) 



dU T (hg) dU(h) 



h=0 



-U f T A . 



(2.83) 



h=0 



We have used the local coordinate version of the right invariant 1-forms and the orthogo- 
nality property of the Ta matrix basis (tr (T^Tg) = — Sab)- 

For simplicity we will assume the following Hamiltonian, when deriving the equations 
of motion: H = \j^ a pl + S(g). Using this Hamiltonian and Eq. 12.821 the equations of 
motion of Eq. |2.78l can be transformed into the simple, well-known form: 



. dU 

U = 5~ < 



dS 



n = -rtr(^dll(r 8 )) = -Ttr{TU 



dH 

as 

dU f) 



dU{r a ) 



-A(U 



nu, 
as 

dU f) 



(2.84) 



with A traceless, antihermitian matrix projector. 

Finally let us calculate in the matrix representation the Jacobian of / phase space func- 
tion, which maps the {U(g),YL(p)} variables to {U(G),Yl(P)}. Lets take the first element 
of the / hypermatrix in Eq. I2.81l and use Eq. I2.82l and Eq. 12 .831 to eliminate the components 
of the right invariant fields: 



Q A (f*r l 



(g) = qHg] 



dG c dUus g / ,air 



dU aB dU^s 
-(U'Ta] 



3a aTT 



7<5 



dgd 

y B u) lS 



(2.85) 



Similar calculation yields the other matrix elements of /: 



Q A (f* 



dp 1 



[U'T A ] 



air a/3 
3a arT 



7<5 



Tb)-y6, 



dP A (f*r l 



an. 



dP A (f* 



an a/3 



y B u) 



yd' 



dp 1 



t a)^^^(T b ) 



7<5- 



(2.86) 
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Chapter 3 

Dynamical staggered f ermions 



It is known for a long time, that for high enough temperatures and/ or densities the quarks 
and gluons are liberated from confinement, the chiral symmetry is restored: the so called 
quark-gluon plasma phase of the matter is created. 

There is a huge literature of this transition: theoretical works based on the symmetries 
of QCD, analytical and numerical calculations in QCD like models and lattice QCD. It is 
worth emphasizing, that the only known way to obtain the properties of the quark-gluon 
plasma from first principles of the theory is lattice QCD. However until recently lattice 
result were usually burdened by large systematical errors: extrapolation to the physical 
quark mass, finite volume effects, missing continuum extrapolations. 

Lattice QCD recently has entered a new era, where we are facing a huge reduction of 
these systematics. In this chapter we describe the details and results of a large scale simu- 
lation, where we attempted to eliminate (almost) all systematics of previous lattice calcu- 
lations. Thus these results can be considered as the final ones, where the only remaining is 
to crosscheck against the work of other groups or against different lattice discretizations. 
The amount of computer work is tremendous, we used O(10 19 ) floating point operations 
on the fastest supercomputers of the world. The algorithmic and theoretical improvements 
still continue, so as the increase in the speed of the computers. We hope that one day these 
results will be just as easy to obtain as getting the value of eg. sin(l.O) using a pocket 
calculator today. 

We emphasize that extensive experimental work is currently being done with heavy 
ion collisions to study the QCD transition (most recently at the Relativistic Heavy Ion 
Collider, RHIC). Moreover there is rich perspective for the future: the heavy ion program is 
expected to start in 2009 at the Large Hadron Collider (LHC) in Geneva and in 2011 at the 
Facility for Ion and Antiproton Research (FAIR) in Darmstadt. Both for the cosmological 
transition and for RHIC, the net baryon densities are quite small, and so the baryonic 
chemical potentials (pi) are much less than the typical hadron masses (~45 MeV at RHIC 
and negligible in the early Universe). A calculation at fi=0 is directly applicable for the 
cosmological transition and most probably also relevant for the transition at RHIC. 
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Let us remark here, that even if the finite temperature equilibrium state of QCD is soon 
going to be solved, there are still important areas with only moderate or no progress. Most 
notably there is the equilibrium state at finite ji, the well-known sign problem prohibited 
calculations for many years. The breakthrough of ||5TI|52|| has opened new possibilities (for 
a recent review of this subfield see 1153 II ), still many questions remain unanswered. 

In Sec. 1 we present the definition of our lattice action, the numerical details of the 
algorithm used for the simulations and finally the concept of line of constant physics (LCP). 
In Sec. 2 we give a detailed list of our simulation points, whereas Sec. 3 is for the physics 
results. 

3.1 Setting up the simulations 

3.1.1 The lattice action 

First we give our definition for the Symanzik improved gauge and for the stout-link im- 
proved fermionic action. We demonstrate that our choice of stout-link improved staggered 
fermionic action has small taste violation, when compared to other staggered actions used 
in the literature to determine the equation of state (EoS) of QCD. 

Isotropic lattice couplings are used, thus the lattice spacings are identical in all direc- 
tions. The lattice action we used has the following form: 

S = S g + S f , (3.1) 

s s = Ef^E^W+^E^fW)' (3-2) 

x 3 }i>v ji^ v 

s f = E{^W^(^ stoM %!/ + ^A,y]- 1/2 ^(y) 

x,y 

+ // s (x)[C<(L[ sto "0x y + m s ^ y ]- 1/ V(y)}, (3.3) 

where W^* , W^* 2 are real parts of the traces of the ordered products of link matrices 
along the 1 x 1, 1 x 2 rectangles in the ji, v plane. The coefficients satisfy Cq + 8c\ = 1 and 
C\ = —1/12 for the tree-level Symanzik improved action. rj u ^ and t] s are the pseudofermion 
fields for u, d and s quarks. ]J(U st0Ut ) is the four-flavor staggered Dirac matrix with stout- 
link improvement |43]|. Let us also note here, that we use the 4th root trick in Eq. (|3.1|) , 
which might lead to problems of locality. 

Our staggered action at a given Nt yields the same limit for the pressure at infinite tem- 
peratures as the standard unimproved action. There are various techniques improving the 
high temperature scaling. However one also has to take into account, that with highly im- 
proved actions (which contain far neighbor interactions) smaller (Nt > 8) lattice spacings 
will not be available. In this case one risks to have large lattice artefacts coming from the 
scale setting procedure. 
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Figure 3.1: Average gauge action densities in various simulations. The blue points are from 
inexact R-algorithm simulations, one can clearly see the stepsize (e) dependence. Previous 
thermodynamical calculations used the R-algorithm, with the stepsize set to the half of 
the light quark mass. The black point is obtained using the RHMC algorithm with large 
stepsize. Note, that it is unsafe to use the R-algorithm for equation of state calculations, 
since the typical size of the subtraction (ie. gauge action difference on zero and non-zero 
temperature lattices) is in the order of the systematical error of the R-algorithm. 



Staggered fermions have an inconvenient property: they violate taste symmetry at fi- 
nite lattice spacing. Among other things this violation results in a splitting in the pion 
spectrum, which should vanish in the continuum limit. The stout-link improvement makes 
the staggered fermion taste symmetry violation small already at moderate lattice spacings. 
We found that a stout-smearing level of N smr =2 and smearing parameter of ,0=0.15 are the 
optimal values of the smearing procedure. 

3.1.2 The algorithm 

The equation of state calculation is an extreme high precision measurement. There are 
many things which can spoil it, one of them is the systematic error coming from the al- 
gorithm. Before our work the R-algorithm was used exclusively for staggered thermody- 
namical calculations, where in principle one has to make an extrapolation in the intrinsic 
parameter of the algorithm (stepsize). These extrapolations were never carried out, in the 
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best cases there were attempts to estimate the systematic errors. The Nt = 6 equation of 
state has forced us to change. Here the measured quantity (action density difference at 
zero and finite temperature lattices) has the same magnitude as the systematic error (see 
Fig. 13.11 ), which is clearly an unsafe situation. Fortunately the development of exact stag- 
gered algorithms were in a prospective phase at that time (rational hybrid Monte-Carlo 
||54| and polynomial hybrid Monte-Carlo |55|). We decided to use the RHMC, the algo- 
rithm which is nowadays obligatory in lattice thermodynamics. It was worth changing, 
the exact RHMC algorithm is significantly faster than the R-algorithm, and one can get rid 
of one systematical error. 

The RHMC technique approximates the fractional powers of the Dirac operator by ra- 
tional functions. Since the condition number of the Dirac operator changes as we change 
the mass, one should determine the optimal rational approximation for each quark mass. 
Note however, that this should be done only once, and the obtained parameters of these 
functions can be used in the entire configuration production. Our choices for the rational 
approximation were as good as few times the machine precision for the whole range of 
the eigenvalues of the Dirac operator. We have also introduced multiple time scales in the 
algorithm ||56l|. 

The time consuming parts of the computations were carried out in single precision. 
This might effect the algorithm in a negative way at two places: firstly the reversibility of 
the trajectory is lost, secondly the precision problems in the accept/ reject step might result 
a bad distribution. 

The reversibility violation is considerably larger than single precision accuracy, even if 
every step was carried out upto single precision through the trajectory. This is due to the 
chaotic nature of the QCD equations of motion. The usual way out is to use double preci- 
sion arithmetics everywhere. However it turns out that at several places single precision 
accuracy is tolerable, if at some critical places high enough precision is chosen. We make 
the force calculation (which is the most time consuming part) in single precision however 
the link and the momentum are updated in a higher precision scheme (in turned out that 
we need at least 80bit precision on larger volumes). In this case the reversibility will be 
exact in single precision. The link and the momentum might differ after going forth and 
back along a trajectory, but only in high precision. The forces will be bit by bit the same 
in the forward and backward directions, which ensures that the two cannot deviate from 
each other. 

For the other problem (leakage of precision in the accept/ reject step) one can use mixed 
precision inverters, which work in single precision for most of the time. Here one adds 
intermediate double precision steps, with which one can achieve even double precision 
accuracy. To be on the safe side on one of our largest lattices we have cross-checked the 
results with a fully double precision calculation, the results were the same within the er- 
rorbars. We were also constantly monitoring the (exp(— AH)) expectation value, and find 
no statistically significant deviation from 1. 

We based our code on the publicly available MILC lattice gauge theory code, how- 



46 



3.1. Setting up the simulations 




Figure 3.2: The line of constant physics. The left panel shows the strange mass as a function 
of j6 along LCP1. LCP1 is an approximate LCP, it was obtained by using the (p and K masses 
(see text). LCP2 is a refinement of LCP1. The right panel shows the strange mass (red) and 
20 times the light quark mass (blue) in lattice units as functions of jS = 6/g 2 along LCP2. 

ever several parts were (re)written by ourselves. Most of the code was written in two 
independent copies, the two versions agree upto machine precision. These (among other 
things) include the staggered matrix multiplication, solvers, smearing and measurement 
routines. We have developed the code for four different architectures (Intel P4, AMD 
Opteron, Nvidia Graphics Card, IBM Blue Gene L), each version required careful opti- 
mization. For some details of the implementations see eg. f57ll3Tf . 

3.1.3 Line of constant physics (LCP) 

Let us discuss the determination of the LCP. The LCP is defined as relationships between 
the bare lattice parameters (/3 and lattice bare quark masses m u ^ and m s ). These relation- 
ships express that the physics (e.g. mass ratios) remains constant, while changing any of 
the parameters. It is important to emphasize that the LCP is unambiguous (independent of 
the physical quantities, which are used to define the above relationships) only in the con- 
tinuum limit (/3 — > oo). For our lattice spacings fixing some relationships to their physical 
values means that some other relationships will slightly deviate from the physical one. In 
thermodynamics the relevance of LCP comes into play when the temperature is changed 
by /3 parameter. Then adjusting the mass parameters (m u & and m s ) is an important issue, 
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neglecting this in simulations can lead to several % error in the EoS [58]. 

A particularly efficient (however only approximate, see later) way to obtain an LCP is 
by using simulations with three degenerate flavors with lattice quark mass m.a. The leading 
order chiral perturbation theory implies the mass relation for ss mesons. The strange quark 
mass is tuned accordingly as 

m 2 PS /m 2 v \ mq=ms = (2m| ~ m \)l m \> ( 3 - 4 ) 

where mpg and my are the pseudoscalar and vector meson masses in the simulations with 
three degenerate quarks. The light quark mass is calculated using the ratio m u & = m s /25 
obtained by experimental mass input in the chiral perturbation theory. We obtain m s (jS) 
as shown in the left panel of Fig. 13.21 This (approximate) line of constant physics is called 
LCP1 later, the equation of state calculations were carried out along this line. 

Our approach using Eq. (|3.4|) is appropriate if in the Hf=2+1 theory the vector meson 
mass depends only weakly on the light quark masses and the chiral perturbation theory 
for meson masses works upto the strange quark mass. After applying the LCP1 we cross- 
checked the obtained spectrum of the rij=2+l simulations. These simulations showed, 
however, that the hadron mass ratios slightly differ from their physical values on the 5- 
10% level. In order to eliminate all uncertainties related to an unphysical spectrum, we 
determined a new line of constant physics. The new LCP (which is called LCP2 afterwards) 
was defined by fixing m^/ /k an d mK/m n to their experimental values (right panel of Fig. 
I3.2|) . The more precise LCP2 was used for simulations to determine the order of the QCD 
transition, and to measure the transition temperature in physical units. 

We have also carried out rij = 2 + 1 flavor T = simulations on LCP2. Chiral extrap- 
olation to the physical pion mass led to m^l /k and m^/ m n values, which agree with the 
experimental numbers on the 2% level. (Differences resulting from various fitting forms 
and finite volume corrections were included in the systematics.) This is the accuracy of 
LCP2. 

In order to be sure that our results are safe from ambiguous determination of the over- 
all scale, and to prove that we are really in the a 2 scaling region, we carried out a con- 
tinuum extrapolation for three additional quantities which could be similarly good to set 
the scale (we normalized them by fx, for fx determination in staggered QCD see |59]]). 
Fig. 13.31 shows the measured values of mx* I fx, fnl fx and r§fx, at different lattice spac- 
ings and their continuum extrapolation. Our three continuum predictions are in complete 
agreement with the experimental results (note, that tq can not be measured directly in ex- 
periments; in this case the original experimental input is the bottonium spectrum which 
was used by the MILC, HPQCD and UKQCD collaborations to calculate tq on the lattice 
B53|6yl). 

It is important to emphasize that at lattice spacings given by Nf=4 and 6 the overall 
scales determined by fx and r$ are differing by ^20-30%, which is most probably true 
for any other staggered formulation used for thermodynamical calculations. Since the 
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Figure 3.3: Scaling of the mass of the K*(892) meson, the pion decay constant and tq to- 
wards the continuum limit. As a continuum value (filled boxes) we took the average of the 
continuum extrapolations obtained using our 2 and our 3 finest lattice spacings. The dif- 
ference was taken as a systematic uncertainty, which is included in the shown errors. The 
quantities are plotted in units of the kaon decay constant. In case of the upper two panels 
the bands indicate the physical values of the ratios and their experimental uncertainties. 
For tq (lowest panel) in the absence of direct experimental results we compare our value 
with the r f K obtained by the MILC, HPQCD and UKQCD collaborations 15911501. 
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determination of the overall scale has a ^20-30% ambiguity, the value of T c can not be 
determined with the required accuracy. 

3.2 Simulation points 

For our thermodynamical calculations we have used two LCPs: LCP1 and LCP2. The 
LCP1 can be considered as an approximate LCP, this was used for the equation of state 
calculation. The LCP2 was determined using the nt = 2 + 1 simulations carried out along 
LCP1, it can be considered as a refinement of LCP1. We used it to determine the order 
of phase transition and the transition temperature. In this section we list the simulation 
points along the two LCPs. 

3.2.1 Along LCP1 

The determination of the EoS needs quite a few simulation points. Results are needed on 
finite temperature lattices (Nf=4 or 6) and on zero temperature lattices (Nt 4 or 6) at 
several f> values (we used 16 different /3 values for Nf =4 and 14 values for Nt=6). Since our 
goal is to determine the EoS for physical quark masses we have to determine quantities in 
this small physical quark mass limit (we call these /3 dependent bare light quark masses 
m ud (phys)). 

For our finite temperature simulations (Nf=4,6) we used physical quark masses. The 
spatial sizes were always at least 3 times the temporal sizes. For the whole /3 range on 
Nt = 4 we checked that by increasing the N s /Nt ratio from 3 to 4 the results remained the 
same within our statistical uncertainties. 

In the chirally broken phase (our zero temperature simulations, thus lattices for which 
Nt 4 or 6, belong always to this class) chiral perturbation theory can be used to extrap- 
olate by a controlled manner to the physical light quark masses. Therefore for most of our 
simulation points^ we used four pion masses (m n ^250, 320, 380 and 430 MeV), which 
were somewhat larger than the physical one. (To simplify our notation in the rest of this 
section we label these points as 3,5,7 and 9 times tn ud (phys).) It turns out that the chiral 
condensates at all the four points can be fitted by linear function of pion mass squared 
with good x 2 - (Later we will show, that only the chiral condensate is to be extrapolated to 
get the EoS at the physical quark mass.) The volumes were chosen in a way, that for three 
out of these four quark masses the spatial extentions of the lattices were approximately 
equal or larger than four times the correlation lengths of the pion channel. We checked 
for a few /3 values that increasing the spatial and/ or temporal extensions of the lattices re- 
sults in the same expectation values within our statistical uncertainties. (For 3-m ud (phys) 
values the spatial lengths of the lattices were only three times the correlation length of the 

1 In the f> = 3.0. .3.4 range the T = simulations were carried out at m u ^(phys). 
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Table 3.1: Summary of our simulation points along LCP1. For the physical light quark 
masses (we call them m uc i(phys)) 25 times smaller values were taken than for the strange 
mass. T^0 simulations were performed with the above m s and pairs, and at 5 different 
m ui i values: {1,3,5,7, 9}-m U( i(phys). T=0 simulations were performed with the above m s and 
j8 pairs, but at 4 different m u ^ values: {3,5,7,9}-m u a(phys). The total number of trajectories 
divided by 100 are collected in the # columns. The left column shows the N;=4, whereas 
the right column shows the Nt=6 data. (For an explanation of our labeling see the text.) 
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pion channel. However, excluding this point from the extrapolations, the results do not 
change.) 

A detailed list of our simulation points at zero and at non-zero temperature lattices are 
summarized in Table |3~T1 

3.2.2 Along LCP2 

In order to perform the necessary renormalizations of the measured quantities and to fix 
the scale in physical units we carried out T = simulations on our new LCP (c.f . Table l3~2l> . 
Six different /3 values were used. Simulations at T=0 with physical pion masses are quite 
expensive and in our case unnecessary (chiral perturbation theory provides a controlled 
approximation at vanishing temperature). Thus, for each /3 value we used four different 
light quark masses, which resulted in pion masses somewhat larger than the physical one 
(the m n values were approximately 250 MeV, 320 MeV, 380 MeV and 430 MeV), whereas 
the strange quark mass was fixed by the LCP at each /3. The lattice sizes were chosen to 
satisfy the m n N s > 4 condition. However, when calculating the systematic uncertainties of 
meson masses and decay constants, we have taken finite size corrections into account using 
continuum finite volume chiral perturbation theory [6T| (these corrections were around or 
less than 1%). We have simulated between 700 and 3000 RHMC trajectories for each point 
in Table [3721 

The T^0 simulations (c.f. Table l33l> were carried out along our LCP (that is at physical 
strange and light quark masses, which correspond to ra^=498 MeV and m n =l35 MeV) at 
four different sets of lattice spacings (Nt = 4, 6, 8 and 10) and on three different volumes 
(N s I Nt was ranging between 3 and 6). We have observed moderate finite volume effects on 
the smallest volumes for quantities which are supposed to depend strongly on light quark 
masses (e.g. chiral susceptibility). To determine the transition point we used N s /Nt > 4, 
for which we did not observe any finite volume effect. The number of RHMC trajectories 
were between 1500 and 8000 for each parameter set (the integrated autocorrelation time 
was smaller or around 10 for all our runs). 
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Table 3.2: Lattice parameters and sizes of our zero temperature simulations. The strange 
quark mass is varied along the LCP as /3 is changed. The light quark masses, listed at each 
(/3,m s ) values, correspond approximately tom n values of 250 MeV, 320 MeV, 380 MeV and 
430 MeV. 



temporal size (N t ) 


/3 range 


spatial sizes (N s ) 


4 


3.20 - 3.50 


12,16,24 


6 


3.45 - 3.75 


18,24,32 


8 


3.57 - 3.76 


24,32,40 


10 


3.63 - 3.86 


28,40,48 



Table 3.3: Summary of the T^0 simulation points. 
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3.3 Improvement over previous results 

As we have already mentioned in the introduction, there are many lattice results on QCD 
thermodynamics. In this section we highlight the points, where we have made improve- 
ments on previous calculations. 

Physical quark masses 

We decided to use physical values for the quark masses. Owing to the computational 
costs this is a great challenge in lattice QCD. Previous analyses used computationally less 
demanding non-physically large quark masses. 

On the one hand, results with Wilson fermions [62ll63| were obtained with pion masses 
m n > 540 MeV when approaching the thermodynamical limit (since lattice QCD can give 
only dimensionless combinations, it is more precise to say that m n /nip > 0.6, where m p is 
the mass of the rho meson). 

On the other hand, in staggered simulations one can afford considerably smaller quark 
masses. The MILC collaboration ||64l |65l is currently using two light quark masses (0.1 
and 0.2 times m s ), the Bielefeld-Brookhaven-Columbia-RIKEN collaboration is studying 
thermodynamics down to a pion mass of ~ 320 MeV on Nf=4 and 6 lattices ||66|. How- 
ever these numbers should be taken with a grain of salt. Staggered fermions suffer from 
taste violation. Therefore there is a large (usually several hundred MeV), unphysical mass 
splitting between this lightest pion state and the higher lying other pion states. This mass 
splitting results in an unphysical spectrum. The artificial pion mass splitting disappears 
only in the continuum limit. For some choices of the actions the restoration of the proper 
spectrum happens only at very small lattice spacings, whereas for other actions somewhat 
larger lattice spacings are already satisfactory. 

The finite temperature transition is related to the spontaneous breaking of the chiral 
symmetry (which is driven by the pion sector) and the three physical pions have masses 
smaller than the transition temperature, thus the numerical value of T c could be sensitive 
to the unphysical spectrum. Furthermore, the order of the transition depends on the quark 
mass. In three-flavor QCD for vanishing quark masses the transition is of first-order. For 
intermediate masses it is most probably a crossover. For infinitely heavy quark masses the 
transition is again first-order. Therefore the physical quark masses should be used directly. 

It is also important to mention that though at T=0 chiral perturbation theory provides a 
technique to extrapolate to physical m n , unfortunately no such controllable method exists 
around T c . 

Continuum limit 

The second ingredient is to remove the uncertainty associated with the lattice discretiza- 
tion. Discretization errors disappear in the continuum limit; however, they strongly influ- 
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ence the results at non-vanishing lattice spacing. 

For lattice spacings which are smaller than some approximate limiting value the di- 
mensionless ratio of different physical quantities have a specific dependence on the lattice 
spacing (for staggered QCD the continuum value is approached in this region by correc- 
tions proportional to the square of the lattice spacing). For these lattice spacings we use the 
expression: a 2 scaling region. Clearly, results for at least three different lattice spacings are 
needed to decide, whether one is already in this scaling region or not (two points can al- 
ways be fitted by Co+C2« 2 , independently of possible large higher order terms). Only using 
the a 2 dependencies in the scaling region, is it possible to unambiguously define the abso- 
lute scale of the system. Outside the scaling regiorH different quantities lead to different 
overall scales, which lead to ambiguous values for e.g. T c . 

In three-flavour unimproved staggered QCD, using a lattice spacing of about 0.28 fm, 
the first-order and the crossover regions are separated by a pseudoscalar mass of m n /C ~ 
300 MeV. Studying the same three-flavour theory with the same lattice spacing, but with an 
improved p4 action (which has different discretization errors) we obtain m n /C ~ 70 MeV. 
In the first approximation, a pseudoscalar mass of 140 MeV (which corresponds to the 
numerical value of the physical pion mass) would be in the first-order transition region, 
whereas using the second approximation, it would be in the crossover region. The different 
discretisation uncertainties are solely responsible for these qualitatively different results 

In summary the proper approach is to extrapolate to vanishing lattice spacings using 
lattices which are alrady in the scaling regime. We approach the scaling region by using 
four different sets of lattice spacing, which are defined as the transition region on Nt =4,6,8 
and 10 lattices. The results show (not surprisingly) that the coarsest lattice with Nf=4 is not 
in the a 2 scaling region, whereas for the other three a reliable continuum limit extrapolation 
can be carried out. In case of the equation of state we only have two lattice spacings (Nt = 4 
and 6), for the continuum limit we have to wait for results on finer lattices. 

Lattice artefacts for T = and for T — > 00 

For the staggered formulation of quarks the physically almost degenerate pion triplet has 
an unphysical non-degeneracy (so-called taste violation). This mass splitting Am 2 van- 
ishes in the continuum limit as a — >0. Due to our smaller lattice spacing and particularly 
due to our stout-link improved action the splitting Am 2 is much smaller than that of the 
previously or currently used staggered actions in thermodynamics. In order to illustrate 
the advantage of the stout-link action Fig. 13.41 compares the taste violation in different 

2 Note, that outside the scaling region even a seemingly small lattice spacing dependence can lead to an 
incorrect result. An infamous example is the Naik action [67] in the Stefan-Boltzmann limit: Nf=4 and 6 are 
consistent with each other with a few % accuracy, but since they are not in the scaling region they are 20% 
off the continuum value. 
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Figure 3.4: Pion mass splitting Ara^ = — as a function of m\ in units of tq for 
different works on lattice thermodynamics with staggered quarks. The lattice spacings are 
the same as those at the finite temperature transition point. The mass of the Goldstone 
pion is denoted by m n , that of the first non-Goldstone mode is by m N c- The vertical black 
line corresponds to the physical value of m\. The taste violation of our stout-link improved 
action is smaller than that of any other actions used in the literature. 



approaches of the literature, which have been used for staggered thermodynamics. Re- 
sults on the pion mass splitting for p4 improved (used by Bielefeld-Brookhaven-Columbia- 
RIKEN collaboration ||66l), asqtad improved (used by MILC collaboration ||64l [65] ) and 
stout-link improved (this work) staggered fermions are shown. The parameters were cho- 
sen to be the ones used by the different collaborations at the finite temperature transition 
point. 

At infinitely large temperatures improved actions (p4 |7T| or asqtad |[72l [73| action) 
show considerably smaller discretizations errors, than the standard staggered action (used 
by this work). However our choice of action is about an order of magnitude faster than 
e.g. p4, we decided to use this less improved action, with which our CPU resources 
made it possible to study several lattice spacings (Nf=4 and 6 for the equation of state 



56 



3.3. Improvement over previous results 



and Nt = 4, 6, 8 and 10 for determining the order of the phase transition and the transition 
temperature). This turned out to be extremely beneficial, when converting the transition 
temperature into physical units. In particular the T = simulations -which are used to do 
this conversion- have very large lattice artefacts at Nt = 4 and 6 lattice spacings and can 
not be used for controlled continuum extrapolations. The high-temperature improvement 
is not designed to reduce these artefacts. 

Setting the physical scale 

An additional problem appears if we want to give dimensionful predictions with a few 
percent accuracy. As we already emphasized lattice QCD predicts dimensionless com- 
binations of physical observables. For dimensionful predictions one calculates an experi- 
mentally known dimensionful quantity, which is used then to set the overall scale. In many 
analyses the overall scale is related to some quantities which strictly speaking do not even 
exist in full QCD (e.g. the mass of the rho eigenstate and the string tension are not well 
defined due to decay or string breaking). A better, though still not satisfactory possibility 
is to use quantities, which are well defined, but can not be measured directly in exper- 
iments. Such a quantity is the heavy quark-antiquark potential (V), or its characteristic 
distances: the tq or r\ parameters of V ||74| (r 2 d 2 V/ dr 2 =1.65 or 1, for yq or Y\, respectively). 
For these quantities intermediate lattice calculations and/ or approximations are needed to 
connect them to measurements. These calculations are based on bottonium spectroscopy. 
This procedure leads to further, unnecessary systematic uncertainties. 

The ultimate solution is to use quantities, which can be measured directly in experi- 
ments and on the lattice. We use the decay constant of the kaon /x=159.8 MeV, which has 
about 1% measurement error. Detailed additional analyses were done by using the mass 
of the K* (892) meson m^*, the pion decay constant f n and the value of tq, which all show 
that we are in the a 2 scaling regime and our choice of overall scale is unambiguous (see 
subsection 13.1.3) . 

Algorithm 

In previous staggered thermodynamics simulations the inexact R-algorithm was used ex- 
clusively to simulate three quark flavours. This algorithm has an intrinsic parameter, the 
stepsize which, similarly to the lattice spacing, has to be extrapolated to zero. None of 
the previous staggered lattice thermodynamic studies carried out this extrapolation. Us- 
ing the R-algorithm without stepsize extrapolation leads to uncontrolled systematic errors. 
Instead of using the approximate R-algorithm this work uses the exact RHMC-algorithm 
(rational hybrid Monte-Carlo) ||54|. 
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3.4 Order of the QCD transition 

The nature of the QCD transition affects our understanding of the Universe's evolution (see 
Ref. [75] for example). In a strong first-order phase transition the quark-gluon plasma su- 
percools before bubbles of hadron gas are formed. These bubbles grow, collide and merge, 
during which gravitational waves could be produced [76|. Bary on-enriched nuggets could 
remain between the bubbles, contributing to dark matter. The hadronic phase is the initial 
condition for nucleosynthesis, so inhomogeneities in this phase could have a strong effect 
on nucleosynthesis \\77\\. As the first-order phase transition weakens, these effects become 
less pronounced. Our calculations provide strong evidence that the QCD transition is a 
crossover and thus the above scenarios — and many others — are ruled out. 

There are some QCD results and model calculations to determine the order of the tran- 
sition at n=0 and ji^O for different fermionic contents (compare refs 1751 1791 1501 ETl 1871 l83l 
|84H85l|86l|871). Unfortunately, none of these approaches can give an unambiguous answer 
for the order of the transition for physical values of the quark masses. The only known 
systematic technique which could give a final answer is lattice QCD. 

There are several lattice results for the order of the QCD transition (for the two most 
popular lattice fermion formulations see refs 1188 II and 1162 II ), although they have unknown 
systematics. As we have already emphasized in the previous section from the lattice point 
of view there are two important 'ingredients' to eliminate these systematic uncertainties: 
one has to use physical quark masses and carry out a continuum extrapolation. 

Our goal is to identify the nature of the transition for physical quark masses as we ap- 
proach the continuum limit. We will study the finite size scaling of the lattice chiral suscep- 
tibilities ,^(N S , Nt)=d 2 / (dm 2 d )(T /V)- logZ, where m u ^ is the mass of the light u,d quarks 
and N s is the spatial extension. This susceptibility shows a pronounced peak around the 
transition temperature (T c ). For a real phase transition the height of the susceptibility peak 
increases and the width of the peak decreases when we increase the volume. For a first- 
order phase transition the finite size scaling is determined by the geometric dimension, 
the height is proportional to V, and the width is proportional to 1 / V . For a second-order 
transition the singular behaviour is given by some power of V, defined by the critical expo- 
nents. The picture would be completely different for an analytic crossover. There would be 
no singular behaviour and the susceptibility peak does not get sharper when we increase 
the volume; instead, its height and width will be V independent for large volumes. 

Fig. 13.51 shows the susceptibilities for the light quarks for Nf =4 and 6, for which we 
used aspect ratios r = N s /Nt ranging from 3 to 6 and 3 to 5, respectively. A clear signal 
for an analytic crossover for both lattice spacings can be seen. However, these curves do 
not say much about the continuum behaviour of the theory. In principle a phenomenon 
as unfortunate as that in the three-flavour theory could occur |[68~f, in which the reduction 
of the discretization effects changed the nature of the transition for a pseudoscalar mass of 
wl40 MeV. 

Because we are interested in genuine temperature effects we subtract the T=0 suscepti- 
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Figure 3.5: Susceptibilities for the light quarks for Nf=4 (left panel) and for Nf=6 (right 
panel) as a function of 6/g 2 , where g is the gauge coupling (T grows with 6/g 2 ). The largest 
volume is eight times bigger than the smallest one, so a first-order phase transition would 
predict a susceptibility peak that is eight times higher (for a second-order phase transition 
the increase would be somewhat less, but still dramatic). Instead of such a significant 
change we do not observe any volume dependence. Error bars are s.e.m. 




Figure 3.6: Normalized susceptibilities T 4 / (m 2 Ax) for the light quarks for aspect ratios 
r=3 (left panel) r=4 (middle panel) and r=5 (right panel) as functions of the lattice spacing. 
Continuum extrapolations are carried out for all three physical volumes and the results are 
given by the leftmost blue diamonds. Error bars are s.e.m with systematic estimates. 
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Figure 3.7: Continuum extrapolated susceptibilities T 4 / (m 2 Ax) as a function of 1/(T^V). 
For true phase transitions the infinite volume extrapolation should be consistent with zero, 
whereas for an analytic crossover the infinite volume extrapolation gives a non-vanishing 
value. The continuum-extrapolated susceptibilities show no phase-transition-like volume 
dependence, though the volume changes by a factor of five. The V^oo extrapolated value 
is 22(2) which is lie away from zero. For illustration, we fit the expected asymptotic 
behaviour for first-order and 0(4) (second order) phase transitions shown by dotted and 
dashed lines, which results in chance probabilities of 10~ 19 (7 x 10 -13 ), respectively. Error 
bars are s.e.m with systematic estimates. 

bility and study only the difference between T^O and T=0 at different lattice spacings. To 
do it properly, when we approach the continuum limit the renormalization of % has to be 
performed. This leads to m 2 Ax, which we study (for the details see subsection |3.5.1|) . 

To give a continuum result for the order of the transition we carry out a finite size 
scaling analysis of the dimensionless quantity T 4 / (tn 2 Ax) directly in the continuum limit. 
For this study we need the height of the susceptibility peaks in the continuum limit for 
fixed physical volumes. The continuum extrapolations are done using four different lattice 
spacings (Nf=4,6,8 and 10). The volumes at different lattice spacings are fixed in units of 
T Cf and thus VT^=3 3 A 3 and 5 3 were chosen. (In three cases the computer architecture did 
not allow us to take the above ratios directly. In these cases, we used the next possible 
volume and interpolated or extrapolated. The height of the peak depends weakly on the 
volume, so these procedures were always safe.) Altogether we used twelve different lattice 
volumes ranging from 4 • 12 3 to 10 • 48 3 at T > 0. For the T = runs lattice volumes from 
24 • 12 3 up to 56 • 28 3 were used. The number of trajectories were between 1500 and 8000 
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for T > and between 1500 and 3000 for T = 0, respectively. Fig. !3.6l shows the continuum 
extrapolation for the three different physical volumes. The N;=4 results are slightly off but 
the Nf=6,8 and 10 results show a good a 2 ocl/N 2 scaling. 

Having obtained the continuum values for T 4 / (m 2 Ax) at fixed physical volumes, we 
study the finite size scaling of the results. Fig. 13.71 shows our final results. The volume 
dependence strongly suggests that there is no true phase transition but only an analytic 
crossover in QCD. 



3.5 Transition temperature 

There are several results in the literature for T c using both staggered and Wilson fermions 
ll89ll90l[62ll9Tll64ll66| . There is an additional limitation of these, which has not been men- 
tioned before. This problem is related to an implicit assumption about a a real singularity, 
thus ignoring the analytic cross-over feature of the finite temperature QCD transition. 

As we have seen before the QCD transition at non-vanishing temperatures is an ana- 
lytic cross-over. Since there is no singular temperature dependence different definitions 
of the transition point lead to different values. The most famous example for this phe- 
nomenon is the water-vapor transition, for which the transition temperature can be de- 
fined by the peaks of dp/dT (temperature derivative of the density) and c p (heat capacity 
at fixed pressure). For pressures (p) somewhat less than p c = 22.064 MPa the transition is 
of first order, whereas at p = p c the transition is second order. In both cases the singular- 
ity guarantees that both definitions of the transition temperature lead to the same result. 
For p > p c the transition is a rapid cross-over, for which e.g. both dp/dT and Cp show 
pronounced peaks as a function of the temperature, however these peaks are at different 
temperature values. Fig. 13.81 shows the phase diagram based on 1192 II . Analogously, there is 
no unique transition temperature in QCD. 

Our goal is to eliminate all the above limitations and give the full answer. We de- 
termine T c using the sharp changes of the temperature (T) dependence of renormalized 
dimensionless quantities obtained from the chiral condensate ({tptp)), quark number sus- 
ceptibility (tiq) and Polyakov loop (P). We expect that all three quantities result in different 
transition points (similarly to the case of the water, c.f. Fig. 13.8ft . 



3.5.1 Chiral susceptibility 

The chiral susceptibility of the light quarks (x) is defined as 

x » =T v&-^ z = (3 ' 5) 

ua ua 

where / is the free energy density. Since both the bare quark mass and the free energy 
density contain divergences, Xftp nas to be renormalized. 
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Figure 3.8: The phase diagram of water around its critical point (CP). For pressures be- 
low the critical value (p c ) the transition is first order, for p > p c values there is a rapid 
crossover. In the crossover region the transition temperatures defined from different quan- 
tities are not necessarily equal. This can be seen for the temperature derivative of the 
density (dp/dT) and the specific heat (c p ). The bands show the experimental uncertainties 
(see HH). 
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The renormalized quark mass can be written as m# M ^ = Z m ■ m u ^. If we apply a mass 
independent renormalization then we have 

2 3 2 2 3 2 
"-S?T = m K-8^' <3 ' 6) 

The free energy has additive, quadratic divergencies. They can be removed by subtracting 
the free energy at T = (this is the usual renormalization procedure for the free energy or 
pressure), which leads to Therefore, we have the following identity: 

« {T) - ' (T = °» = <-4b /R<T> ' <3 - 7) 

ua K,ua 

the right hand side contains only renormalized quantities, which can be determined by 
measuring the susceptibilities of the left hand side (for the above expression we use the 
shorthand notation m 2 d ■ AXipib)- In order to obtain a dimensionless quantity it is natural to 
normalize the above quantity by T 4 (which minimizes the final errors). Alternatively, one 
can use combinations of T and/ or m n to construct dimensionless quantities (though these 
conventions lead to larger errors). Since the transition is a cross-over (c.f. discussion d 
of our Introduction) the maxima of m^/m^ • Ax^/T 2 or ra^/ra^. • Axff give somewhat 
different values for T c . 

The upper panel of Fig. 13.91 shows the temperature dependence of the renormalized 
chiral susceptibility for different temporal extensions (Nf=4,6,8 and 10). For small enough 
lattice spacings, thus close to the continuum limit, these curves should coincide. As it 
can be seen, the Nt = 4 result has considerable lattice artefacts, however the two smallest 
lattice spacings (Nt = 8 and 10) are already consistent with each other, suggesting that 
they are also consistent with the continuum limit extrapolation (indicated by the orange 
band). The curves exhibit pronounced peaks. We define the transition temperatures by 
the position of these peaks. We fitted a second order expression to the peak to obtain its 
position. The slight change due to the variation of the fitting range is taken as a systematic 
error. The left panel of Fig. 13.101 shows the transition temperatures in physical units for 
different lattice spacings obtained from the chiral susceptibility. As it can be seen Nf=6,8 
and 10 are already in the scaling region, thus a safe continuum extrapolation can be carried 
out. The extrapolations based on Nt = 6, 8, 10 fit and Nt = 8, 10 fit are consistent with each 
other. For our final result we use the average of these two fit results (the difference between 
them are added to our systematic uncertainty). Our T=0 simulations resulted in a 2% error 
on the overall scale. Our final result for the transition temperature based on the chiral 
susceptibility reads: 

T c U # ) = 151(3)(3)MeV 7 (3.8) 
where the first error comes from the T^0, the second from the T=0 analyses. 
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Figure 3.9: Temperature dependence of the renormalized chiral susceptibility 
(m 2 Ax x p t p/T 4: ), the strange quark number susceptibility (Xs/T 2 ) and the renormalized 
Polyakov-loop (Pr) in the transition region. The different symbols show the results for 
N t = 4,6,8 and 10 lattice spacings (filled and empty boxes for Nt = 4 and 6, filled and 
open circles for Nt = 8 and 10). The vertical bands indicate the corresponding transition 
temperatures and its uncertainties coming from the T^O analyses. This error is given by 
the number in the first parenthesis, whereas the error of the overall scale determination is 
indicated by the number in the second parenthesis. The orange bands show our contin- 
uum limit estimates for the three renormalized quantities as a function of the temperature 
with their uncertainties. 
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We use the second derivative of the chiral susceptibility (x") at the peak position to 
estimate the width of the peak ((AT C ) 2 = — x(Tc)/x"(Tc))- For the continuum extrapolated 
width we obtained: 

AT c (x # )=28(5)(l)MeV. (3.9) 

Note, that for a real phase transition (first or second order), the peak would have a 
vanishing width (in the thermodynamic limit), yielding a unique value for the transition 
temperature (which then would be called critical temperature). Due to the crossover na- 
ture of the transition there is no such value, there is a range (151 ± 28 MeV) where the 
transition phenomena takes place. Other quantities than the chiral susceptibility could 
result in transition temperatures within this range. 

The MILC collaboration also reported a continuum result on the transition temperature 
based on the chiral susceptibility ||64| . Their result is 169(12)(4) MeV. Note, that their lattice 
spacings were not as small as ours (they used Nf=4,6 and 8), their aspect ratio was quite 
small (N s /Nt=2), they used non-physical quark masses (their smallest pion mass at T^O 
was ~220 MeV), the non-exact R-algorithm was applied for the simulations and they did 
not use the renormalized susceptibility, but they looked for the peak in the bare Xtpxp/T 2 - 
Using T 4 as a normalization prescription (as we did) the transition temperature would 
decrease their T c values by approximately 9 MeV. Note, that their continuum extrapolation 
resulted in a quite large error. Taking into account their uncertainties our result and their 
result agree on the 1-sigma level. 



3.5.2 Quark number susceptibility 

For heavy-ion experiments the quark number susceptibilities are quite useful, since they 
could be related to event-by-event fluctuations. Our second transition temperature is ob- 
tained from the strange quark number susceptibility, which is defined via [|64| 

J_ 9 2 logZ 
T 2 TV dji 2 

where }i s is the strange quark chemical potential (in lattice units). Quark number suscep- 
tibilities have the convenient property, that they automatically have a proper continuum 
limit, there is no need for renormalization. 

The middle panel of Fig. 13.91 shows the temperature dependence of the strange quark 
number susceptibility for different temporal extensions (Nf=4,6,8 and 10). For small enough 
lattice spacings, thus close to the continuum limit, these curves should coincide again (our 
continuum limit estimate is indicated by the orange band). 

As it can be seen, the Nt = 4 results are quite off, however the two smallest lattice 
spacings (Nt = 8 and 10) are already consistent with each other, suggesting that they are 
also consistent with the continuum limit extrapolation. This feature indicates, that they are 
closer to the continuum result than our statistical uncertainty. 



(3.10) 
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Figure 3.10: Continuum limit of the transition temperatures obtained from the renormal- 
ized chiral susceptibility (m^Axf^/l 4 ), strange quark number susceptibility (Xs/T 2 ) and 
renormalized Polyakov-loop (Pr). 



We defined the transition temperature as the peak in the temperature derivative of the 
strange quark number susceptibility, that is the inflection point of the susceptibility curve. 
The position was determined by two independent ways, which yielded the same result. In 
the first case we fitted a cubic polynomial on the susceptibility curve, while in the second 
case we determined the temperature derivative numerically from neighboring points and 
fitted a quadratic expression to the peak. The slight change due to the variation of the 
fitting range is taken as a systematic error. The middle panel of Fig. 13.101 shows the transi- 
tion temperatures in physical units for different lattice spacings obtained from the strange 
quark number susceptibility. As it can be seen Nf=6,8 and 10 are already in the a 2 scaling 
region, thus a safe continuum extrapolation can be carried out. The extrapolations based 
on Nt = 6, 8, 10 fit and Nt = 8, 10 fit are consistent with each other. For our final result 
we use the average of these two fit results (the difference between them is added to our 
systematic uncertainty). The continuum extrapolated value for the transition temperature 
based on the strange quark number susceptibility is significantly higher than the one from 
the chiral susceptibility. The difference is 24(4) MeV. For the transition temperature in the 
continuum limit one gets: 

T c ( Xs ) =175(2) (4) MeV, (3.11) 

where the first (second) error is from the T^O (T=0) temperature analysis (note, that due 
to the uncertainty of the overall scale, the difference is more precisely determined than the 
uncertainties of T c (Xf\p) an d T c (xs) would suggest). El Similarly to the chiral susceptibility 

3 A continuum extrapolation using only the two coarsest lattices (Nt = 4 and 6) yielded T c ~ 190 MeV 
l93ll , where an approximate LCP (LCP1) was used, if the lattice spacing is set by Tq. 
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analysis, the curvature at the peak can be used to define a width for the transition. 

AT c ( Xs ) =42(4)(l)MeV. (3.12) 

3.5.3 Polyakov loop 

In pure gauge theory the order parameter of the confinement transition is the Polyakov- 
loop: 

P = ^3 E tr [ W 4(x,0)*i 4 (x, 1) . . • <i 4 (x,Nf - 1)]. (3.13) 

P acquires a non-vanishing expectation value in the deconfined phase, signaling the spon- 
taneous breakdown of the Z(3) symmetry. When fermions are present in the system, the 
physical interpretation of the Polyakov-loop expectation value is more complicated (see 
e.g.. [94J). However, its absolute value can be related to the quark-antiquark free energy at 
infinite separation: 

|(P)| 2 = exp(-AF w -(r oo)/T). (3.14) 

AFqcj is the difference of the free energies of the quark-gluon plasma with and without the 
quark-antiquark pair. 

The absolute value of the Polyakov-loop vanishes in the continuum limit. It needs 
renormalization. This can be done by renormalizing the free energy of the quark-antiquark 
pair 1195 II . Note, that QCD at T^O has only the ultraviolet divergencies which are already 
present at T=0. In order to remove these divergencies at a given lattice spacing we used a 
simple renormalization condition 1196 II : 

V R (r )=0, (3.15) 

where the potential is measured at T=0 from Wilson-loops. The above condition fixes the 
additive term in the potential at a given lattice spacing. This additive term can be used at 
the same lattice spacings for the potential obtained from Polyakov loops, or equivalently 
it can be built in into the definition of the renormalized Polyakov-loop. 

\(P R )\ = |(P)|exp(y(r )/(2T)), (3.16) 

where V(tq) is the unrenormalized potential obtained from Wilson-loops. 

The lower panel of Fig. 13.91 shows the temperature dependence of the renormalized 
Polyakov-loops for different temporal extensions (Nf=4,6,8 and 10). The two smallest lat- 
tice spacings (Nt = 8 and 10) are approximately in 1-sigma agreement (our continuum 
limit estimate is indicated by the orange band). 

Similarly to the strange quark susceptibility case we defined the transition temperature 
as the peak in the temperature derivative of the Polyakov-loop, that is the inflection point 
of the Polyakov-loop curve. To locate this point and determine its uncertainties we used 
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the same two methods, which were used to determine T c (xs)- The right panel of Fig. 13.101 
shows the transition temperatures in physical units for different lattice spacings obtained 
from the Polyakov-loop. As it can be seen Nt=6,8 and 10 are already in the scaling region, 
thus a safe continuum extrapolation can be carried out. The extrapolation and the deter- 
mination of the systematic error were done as for T c (x s ). The continuum extrapolated 
value for the transition temperature based on the renormalized Polyakov-loop is signifi- 
cantly higher than the one from the chiral susceptibility. The difference is 25(4) MeV. For 
the transition temperature in the continuum limit one gets: 



where the first (second) error is from the T^O (T=0) temperature analysis (again, due to 
the uncertainties of the overall scale, the difference is more precisely determined than the 
uncertainties of T c {x) an d T C (P) suggest). Similarly to the chiral susceptibility analysis, 
the curvature at the peak can be used to define a width for the transition. 



3.5.4 Comparison with the recent result of the BBCR collaboration 

Let us comment here on an independent study on T c based on large scale simulations of the 
Bielefeld-Brookhaven-Columbia-RIKEN group II66II . The p4fat3 action was used, which is 
designed to give very good results in the (T^oo) Stefan-Boltzmann limit (their action is 
not optimized at T=0, which is needed e.g. to set the scale). The overall scale was set 
by tq. The T c analysis based on the chiral susceptibility peak gave in the continuum limit 
T c (^)=192(7)(4) MeV. (The second error, 4 MeV, estimates the uncertainty of the continuum 
limit extrapolation, which we do not use in the following, since we attempt to give a more 
reliable estimate on that.) This result is in obvious contradiction with our continuum result 
from the same observable, which is T c (^)=151(3)(3) MeV. For the same quantity (position 
of chiral susceptibility peak with physical quark masses in the continuum limit) one should 
obtain the same numerical result independently of the lattice action. Since the chance 
probability that we are faced with a statistical fluctuation and both of the results are correct 
is small, we attempted to understand the origin of the discrepancy. We repeated some 
of their simulations and analyses. In these cases a complete agreement was found. In 
addition to their T=0 analyses we carried out an fa determination, too. This /x was used 
to extend their work, to use an LCP based on fx and to determine T c in physical units. 

We summarize the origin of the contradiction between our findings and theirs. The ma- 
jor part of the difference can be explained by the fact, that the lattice spacings of 1166 II are too 
large (>0.20 fm), thus they are not in the a 2 scaling regime, in which a justified continuum 
extrapolation could have been done. Setting the scale by different dimensionful quantities 
should lead to the same result. However at their lattice spacings the overall scales ob- 
tained by tq or by /x can differ by >20%, and even the continuum extrapolated ro/x value 



T C {P) =176 (3) (4) MeV, 



(3.17) 



AT C (P) = 38(5) (1) MeV. 



(3.18) 
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Figure 3.11: Resolving the discrepancy between the transition temperature of Ref. Il66l and 
that of the present work (see text). The major part of the difference can be traced back to 
the unreliable continuum extrapolation of II66II . Left panel: In Ref. II66II tq was used for 
scale setting (filled boxes), however using the kaon decay constant (empty boxes) leads 
to different transition temperatures even after performing the continuum extrapolation. 
Right panel: in our work the extrapolations based on the finer lattices are safe, using the 
two different scale setting methods one obtains consistent results. 



of these scales is about 4-5(7 away from the value given by the literature Il59ll60l . This scale 
ambiguity appears in T c , too (though other uncertainties of II66II , e.g. coming from the de- 
termination of the peak-position, somewhat hide its high statistical significance). We used 
their T c values fixed by their tq scale, and in addition we converted their peak position of 
the chiral susceptibility to T c setting the scale by (In order to ensure the possibility of a 
consistent continuum limit -independently of the actual physical value of Tq- we used for 
both tq and fx the results of B591I601 as Ref. 1166 II did it for tq.) Setting the overall scale by 
predicts a much smaller T c at their lattice spacings than doing it by tq (see left panel of Fig. 
I3.11|) . Even after carrying out the continuum extrapolation the difference does not vanish 
(~ 30 MeV), which means that the lattice spacings >0.20 fm used by ||66| are not in the 
scaling regime. Thus, results obtained with their lattice spacings can not give a consistent 
continuum limit for T c . 

In our case not only Nt = 4 and 6 temporal extensions were used, but more realistic 
Nt=8 and 10 simulations were carried out, which led to smaller lattice spacings. These 
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calculations are already in the a scaling regime and a safe continuum extrapolation can be 
done. For our lattice spacing different scale setting methods give consistent results. This 
is shown on the right panel of Fig. 13.111 (independently of the scale setting one obtains 
the same T c ) and also justified with high accuracy by Fig. 13.31 where tq/k converges to the 
physical value on our finer lattices. As it can be seen on the plot, using only our Nt = 4 and 
6 results would also give an inconsistent continuum limit. This emphasizes our conclusion 
that lattice spacings > 0.20 fm can not be used for consistent continuum extrapolations. 

The second, minor part of the difference comes from the different definitions of the 
transition temperatures related to the chiral susceptibility. We use the renormalized chiral 
susceptibility with T 4 normalization to obtain the peak position, which yields ~ 9 MeV 
smaller transition temperature than the bare susceptibility normalized by T 2 of Ref. II66II . 



3.6 Equation of state 

The equilibrium description as a function of the temperature is given by the equation of 
state (EoS). The complete determination of the EoS needs non-perturbative inputs, out of 
which the lattice simulation is the most systematic approach. 

The EoS has been determined in the continuum limit for the pure gauge theory 11971 l98l 
199 II . In this case -quenched simulations- the simulations are particularly easy as there is 
no fermionic degree of freedom. The simulated systems is equivalent to that where all the 
fermions are infinitely heavy, thus, is far from the physical situation. (Note, however, that 
even in this relatively simple case there is still a few % difference between the different 
approaches.) 

The situation in the unquenched case (QCD with dynamical quarks) is a bit barouqe. 
There are many results with different flavor content, fermion formulation, quark masses. 
None of them have used the proper physical quark content and none of them has at- 
tempted to carry out a continuum extrapolation. Moreover staggered fermion studies were 
always using inexact R-algorithm (which can result in uncontrolled systematics, see sub- 
section [3X2). Let us give here a brief review of the literature. 



There are published results for two-flavor QCD using unimproved staggered |100l 
I101L and improved Wilson fermions II63L 



There are results available for the 2+1 flavour case, among which the study done by 
Karsch, Laermann and Peikert in the year 2000 ||102 | | , using p4-improved staggered 



fermions, has often been used as the best result of EoS from lattice QCD. An ad- 
ditional drawback of this result, that the concept of the LCP was ignored from the 
calculation. If Karsch,Laermann and Peikert had cooled down e.g. two of their sys- 
tems one at T=3T C and one at T=0.7T C down to T=0, the first system would have had 
approximately 4 times larger quark masses -two times larger pion masses- than the 



70 



3.6. Equation of state 



second one; this unphysical choice is known to lead systematics, which are compa- 
rable to the difference between the interacting and non-interacting plasma 

• Recently the MILC collaboration studied the equation of state along LCP's with two 
light quark masses (0.1 and 0.2 times m s ) at Nj=4 and 6 lattices using asqtad improved 
staggered fermionic action II65II . 

There are ongoing thermodynamics projects improving on previous results. The joint 
Bielefeld-Brookhaven-Columbia-MILC-RIKEN collaboration (hotQCD collaboration) have 
presented new results at the lattice conference | |103| [. Still the results are only available for 



two different lattice spacings (Nt = 4 and 6) and for unphyiscal quark masses. 

There are two important further issues with the equation of state, which have been 
overlooked in previous calculations. Firstly the heavier quarks can have significant con- 
tribution to the pressure even at few times the transition temperature (exploratory inves- 
tigations 11041 show a significant jump in the pressure around ~ 2 • T c due to the charm 
quark). Secondly there is no serious practical obstacle to extend the EoS calculations well 
beyond the usual 4 — 5 • T c 111051 . This opens the possibility to find the missing connection 
between conventional perturbation theory and nonperturbative methods. 

In this section we present our first step towards the final solution of the EoS with proper 
physical quark content. Though we have the EoS on two different sets of lattice spacings 
(Nt = 4 and 6) and one might attempt to do a continuum extrapolation, it is fair to say that 
another set of lattice spacings is needed (Nt=8). One of the reasons is, that in the hadronic 
phase, where the integration for the pressure starts, the lattice spacing is larger than 0.3 fm. 
In this region the lattice artefacts can not be really controlled (and in this deeply hadronic 
case it does not really help that an action is very good at asymptotically high temperatures 
in the free non-interacting gas limit). 



3.6.1 Integral technique 

We shortly review the integral technique to obtain the pressure ||106]| . For large homoge- 
neous systems the pressure is proportional to the logarithm of the partition function: 

pa * = vh l ° sZ{T ' V) = ^ lo SZ(N s ,N f;i 6,m ? ). (3.19) 

' f S 

(Index 'q' refers to the ud and s flavors.) The volume and temperature are connected to the 
spatial and temporal extensions of the lattice: 

V = (N s a)\ T= ^' (3 " 20) 

The divergent zero-point energy has to be removed by subtracting the zero temperature 
(Nt — > oo) part of Eq. (|3.19|) . In practice the zero temperature subtraction is performed by 
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using lattices with finite, but large Nt (called Nto, see Table 13 So the normalized pressure 
becomes: 



T 4 - JV f 



1 1 

j^\ogZ(N s ,Nf,p,m q ) - ——logZiNgOtNtoifrm,,, 



(3.21) 



With usual Monte-Carlo techniques one cannot measure log Z directly, but only its deriva- 
tives with respect to the bare parameters of the lattice action. Having determined the 
partial derivatives one integrates in the multi-dimensional parameter space: 



v_ 

rp4 



Nf 



d{fi, rrin 



(j8o/W?o) 



1 



N t Nl 



dlogZ/dfi 
d\ogZ/dm q 



1 



aiogZo/3/S 
N f0 N s 3 V d\o§Z /dm q 



, (3.22) 



where Z/Zq are shorthand notations for Z(N s ,Nt) /Z(N s o,Nto)- Since the integrand is 
a gradient, the result is by definition independent of the integration path. We need the 
pressure along the LCP, thus it is convenient to measure the derivatives of logZ along 
the LCP and perform the integration over this line in the /3, m ud and m s parameter space. 
The lower limits of the integrations (indicated by /?o and ra^o) were set sufficiently below 
the transition point. By this choice the pressure gets independent of the starting point (in 
other words it vanishes at small temperatures). In the case of 2 + 1 flavor staggered QCD 
the derivatives of log Z with respect to /3 and m q are proportional to the expectation value 
of the gauge action ({S g ) c.f. Eq. (|3.1[> ) and to the chiral condensates ({tfipq)), respectively. 
Eq. p.22|) can be rewritten appropriately and the pressure is given by (in this formula we 
write out explicitely the flavors): 



^ = Nf 



(P,m ud ,m s ) 



d(p,m ud/ m s 



{Po,m udo ,m s0 ) 



1 



N t Nf 



;-s g /fi) 

(#ud) 
(#s) 



1 



:-s g /fi)o 

(#ud)o 
(#s)o 



(3.23) 

where (• • • )o means averaging on a N s 3 • N t o lattice. 

The integral method was originally introduced for the pure gauge case, for which the 
integral is one dimensional, it is performed along the /3 axis. Many previous studies for 
staggered dynamical QCD (e.g. |101l 11071 1102| ) used a one-dimensional parameter space 
instead of performing it along the LCP. Note, that for full QCD the integration should be 
performed along a LCP path in a multi-dimensional parameter space. 

Using appropriate thermodynamical relations one can obtain any thermal properties 
of the system. For example the energy density (e), entropy density (s) and speed of sound 
(Cg) can be derived as 



e = T{dp/dT) - p, 



s = (e + p)T, 



dp 
de 



-j-. (3.24) 
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N t 


p/T 4 


c 1 


X/T 2 


4 


9.12 


1/3 


2.24 


6 


7.86 


1/3 


1.86 


00 


5.21 


1/3 


1 



Table 3.4: Summary of the results for the 2+1 flavor pressure, speed of sound and 1 fla- 
vor quark number susceptibility in the non-interacting Stefan-Boltzmann limit. e/T 4 is 3 
times, whereas s/T 3 is 4 times the normalized value of the pressure (p/T 4 ) in the Stefan- 
Boltzmann limit. The first two lines gives the results for Nf=4,6 and the third line contains 
the results in the continuum (in the thermodynamic limit). 

To be able to do theses derivatives one has to know the temperature along the LCP. Since 
the temperature is connected to the lattice spacing as T = (Nffl) -1 , we need a reliable 
estimate on a. The lattice spacings at different points of the LCP are determined by first 
matching the static potentials for different /3 values at an intermediate distance for m u ^ = 
{3,5}m uc i(phys) quark masses, then extrapolating the results to the physical quark mass. 
Relating these distances to physical observables (determining the overall scale in physical 
units) will be the topic of a subsequent publication. We show the results as a function 
of T/T c . The transition temperature (T c ) is defined by the inflection point of the isospin 
number susceptibility (xi, see later). 

To get the energy density the literature usually uses another quantity, namely e-3p, 
which can be also directly measured on the lattice. In our analysis it turned out to be more 
appropriate to calculate first the pressure directly from the raw lattice data (Eq. (|3.23|> ) and 
then determine the energy density and other quantities from the pressure (Eq. (|3.24|> ). The 
reasons for that can be summarized as follows. As we discussed we perform T^O simula- 
tions with physical quark masses, whereas the subtraction terms from T=0 simulations are 
extrapolated from larger quark masses. This sort of extrapolation is adequate for the chiral 
condensates, for which chiral perturbation techniques work well. Thus, one can choose an 
integration path for the T=0 part of the pressure, which moves along a LCP at some larger 
m ud ( e -g- 9 times m u d{phys)) and then at fixed /3 goes down to the physical quark mass. No 
comparable analogous technique is available for the combination e-3p. 

We have also calculated the pressure for the larger quark masses. Plotting it as a func- 
tion of the temperature the differences between them are significant. As a function oiT/T c 
these differences are smaller, but still remain statistically significant in the 1.2...2.0T C region. 
Note that statements on the mass dependence are only qualitative since such an analysis 
requires the careful matching of the scales at different quark masses. 
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Table 3.5: Numerical values of the pressure for all of our simulation points. The left 
column shows the N;=4, whereas the right column shows the Nt=6 data. Both the raw 
values and the ones scaled by c con t/ 'cn, are given. 



3.6.2 Physics results 

Let us present the results. In order to show how the different quantities scale with the 
lattice spacing we show always N^=4,6 results on the same plot. In addition, in order to 
make the relationship with the continuum limit more transparent we multiply the raw 
lattice results at finite temporal extensions (Nf=4,6) with c CO nt/c^ t , where the c values are 
the results in the free non-interacting plasma (Stefan-Boltzmann limit). These c values 
are summarized in Table 13.6.11 for the pressure, speed of sound, and for the quark num- 
ber susceptibility at Nf=4,6 and in the continuum limit. By this multiplication the lattice 
thermodynamic quantities should approach the continuum Stefan-Boltzmann values for 
extreme large temperatures. 

Table 13.51 contains our most important numerical results. We tabulated the raw and 
normalized pressure values for both lattice spacings and for all of our simulation points. 
This data set and Eq. (|3.24[) were used to obtain the following figures. Fig. 13.121 shows 
the equation of state on Nt=4,6 lattices. The pressure (left panel) and e (right panel) are 
presented as a function of the temperature. The Stefan-Boltzmann limit is also shown. Fig. 
!3.13l shows the entropy density (left panel) and the speed of sound (right panel), which can 
be obtained by using the pressure and energy density data (c.f. sT=e+p and c^=dp/de) of 
the previous Fig. 13.121 Clearly, the uncertainties of the pressure and those of the energy 
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Figure 3.12: a.) The left panel shows the pressure p, as a function of the temperature. Both 
Nt=4: (red, upper curve) and Nt=6 (blue, lower curve) data are obtained along the LCP. 
They are normalized by T 4 and scaled by c CO nt/cN t (see text and Table 13371) . In order to 
lead the eye lines connect the data points, b.) The right panel is the energy density (e), 
red (upper) and blue (lower) for Nf =4 and 6 respectively. This result was obtained directly 
from the pressure. 



density cumulate in the speed of sound, therefore it is less precisely determined. 



Light and strange quark number susceptibilities (x u d an d Xs) are defined via ||64 | 

Xq _ N t 3 2 logZ 



T 2 



N s 3 dp 2 



(3.25) 



}i q =0 



where p u ^ and p s are the light and strange quark chemical potentials (in lattice units). With 
the help of the quark number operators 



i a 



the susceptibilities can be written as 



logdet($ + m,a), 



Xq 

T2 



Nt_ 

N c 3 



{Qq)jl q =0 + 



The first term is usually referred as disconnected, the second as connected part. The con- 
nected part of the light quark number susceptibility is 2 times the susceptibility of the 
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Figure 3.13: The entropy density (upper left panel, normalized by T 3 ) and the speed of 
sound (upper right panel). The isospin susceptibility (lower left panel, normalized by T 2 ) 
and the connected part of the strangeness susceptibility (lower right panel, normalized by 
T 2 ). The labeling is the same as for Fig. 13.121 



isospin number (xi). It is presented on the left panel of Fig. 13.131 For our statistics and 
evaluation method the disconnected parts are all consistent with zero and their value is 
far smaller than those of the connected parts. The right panel of Fig. 13.131 contains the 
connected part of the strange number susceptibility. 
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Chapter 4 
Summary 



The focus of this work was implementing and applying dynamical fermions in lattice 
QCD. 

The first part was exploring the unknown territory of dynamical algorithms for the 
overlap fermion. The conventional dynamical algorithm (Hybrid Monte Carlo) fails to 
work for this type of fermion. The failure was identified, it is due to the change in the 
topological charge. We have proposed a possible workaround for this problem, and shown 
that with this modification the algorithm works reasonably. There were further modifica- 
tions necessary to increase the performance to an acceptable level. At the end we have 
determined the topological susceptibility as the function of quark mass. The result was the 
first in lattice QCD which has shown the suppression of the susceptibility for small quark 
masses at finite lattice spacing. The simulations however considerably more expensive 
than the ones with other fermion formulations. Therefore it still seems to be more benefi- 
cial to investigate the current algorithms or to come up with new ones than starting larger 
scale physics projects with dynamical overlap fermions. 

The second part of the work was the determination of bulk properties of the finite tem- 
perature QCD matter using dynamical improved staggered fermions. This was a large 
scale project and attempted to give final answers based on a first principles approach. In 
order to achieve our goal we have done the simulations for physical values of the quark 
masses and carried out a continuum extrapolation wherever possible. Firstly we have 
determined the nature of the transition based on a finite size scaling analysis of a suscep- 
tibility type quantity. The transition turned out to be a smooth crossover, that is no sign of 
singularity has been found in the thermodynamical limit. Secondly we have determined 
some typical temperatures in physical units where this crossover takes place. Since there 
is no singularity, there is no unique temperature which can be identified as a critical tem- 
perature. We have calculated transition temperatures from various quantities: peak of 
the chiral susceptibility, inflection point of the quark number susceptibility and inflection 
point of the Polyakov loop. The second two gave significantly higher temperature val- 
ues, than the first one. Thirdly we have presented result on the equation of state towards 
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the continuum limit. From the currently available lattice spacings a reliable continuum 
extrapolation cannot be carried out, this is left for the future. 
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